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ABSTRACT 


The  problem  of  coupled  roll,  sway,  and  yaw  stability  analysis  of  submersible 
vehicles  is  analyzed,  with  particular  emphasis  on  nonlinear  studies.  Previous 
results  had  indicated  that  a  primary  loss  of  stability  is  through  the  devel¬ 
opment  of  limit  cycles.  This  loss  of  stability  is  due  to  the  coupling  of  roll 
into  sway  and  yaw  and  cannot  be  predicted  by  considering  the  uncoupled  dy¬ 
namics.  In  this  study,  it  is  shown  that  the  mechanism  of  loss  of  stability  is 
through  bifurcations  to  periodic  solutions.  These  are  characterized  as  either 
subcritical  or  supercritical,  depending  on  the  sign  of  a  certain  nonlinear  coeffi¬ 
cient.  Implications  of  these  results  to  vehicle  performance  and  operations  are 
discussed. 
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I.  INTRODUCTION 


A.  PROBLEM  STATEMENT 

The  dynamic  response  of  a  submersible  vehicle  operating  at  the  extremes 
of  its  operational  envelope  is  becoming  increasingly  important  in  order  to  en¬ 
hance  vehicle  operations.  Traditionally,  dynamic  stability  of  motion  is  studied 
using  eigenvalue  analysis  where  the  equations  of  motion  are  linearized  around 
nominal  straight  line  level  flight  paths  [Arentzen  &  Mandel  (I960)],  [Clayton 
&  Bishop  (1982)],  [Feldman  (1987)].  Directional  stability  in  the  horizontal 
plane  is  normally  studied  assuming  that  coupling  between  sway/yaw  and  roll 
does  not  exist.  Relaxing  this  approximation  can  lead  to  an  oscillatory  loss  of 
directional  stability  [Cunningham  (1993)]  which  cannot  be  predicted  by  un¬ 
coupled  sway/yaw  motions  only.  This  oscillatory  loss  of  stability  can  generate 
limit  cycles  in  the  system,  as  was  confirmed  numerically  in  previous  studies 
[Cunningham  (1993)].  In  order  to  gain  a  better  understanding  of  the  mecha¬ 
nism  of  this  type  of  loss  of  stability  and  the  stability  properties  of  the  resulting 
limit  cycles,  it  is  necessary  to  perform  a  systematic  nonlinear  analysis  which 
is  precisely  the  scope  of  this  work. 
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B.  OBJECTIVES  AND  OUTLINE 


In  this  work  we  examine  the  problem  of  stability  of  motion  with  controls 
fixed  in  the  horizontal  plane,  with  particular  emphasis  on  the  mechanism  of 
loss  of  stability  of  straight  line  motion.  Coupling  between  sway /yaw  and  roll 
is  taken  into  consideration.  This  has  its  origins  in  both  inertial  and  hydro- 
dynamic  coupling.  We  concentrate  on  an  oscillatory  loss  of  stability  case, 
where  one  pair  of  complex  conjugate  of  eigenvalues  of  the  linearized  system 
matrix  crosses  the  imaginary  axis.  This  loss  of  stability  occurs  in  the  form 
of  generic  bifurcations  to  periodic  solutions  [Guckenheimer  &  Holmes  (1983)]. 
Taylor  expansions  and  center  manifold  approximations  are  employed  in  or¬ 
der  to  isolate  the  main  nonlinear  terms  that  influence  system  response  after 
the  initial  loss  of  stability  [Hassard  &  Wan  (1978)].  Integral  averaging  is 
performed  in  order  to  combine  the  nonlinear  terms  into  a  design  stability  co¬ 
efficient  [Chow  &  Mallet-Paret  (1977)].  Special  attention  is  paid  to  the  study 
of  the  quadratic  drag  terms  as  they  constitute  some  of  the  main  nonlinear 
terms  of  the  equations  of  motion.  The  difficulty  associated  with  the  nons¬ 
moothness  of  the  absolute  value  nonlinearities  is  dealt  with  by  employing  the 
concept  of  generalized  gradient  [Clarke  (1983)],  a  technique  which  was  utilized 
in  [Papadimitriou  (1994)].  This  has  the  advantage  of  keeping  the  hnear  terms 
constant,  unlike  the  linear/cubic  approximation  typically  used  in  ship  roll  mo¬ 
tion  studies  [Dalzell  (1978)],  where  the  linear  damping  coefficient  is  a  function 
of  the  assumed  amplitude  of  motion. 

Vehicle  modeling  in  this  work  follows  standard  notation  [Gertler  &  Hagen 
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(1976)],  [Smith  et  al  (1978)],  and  numerical  results  are  presented  for  a  variant 
of  the  Swimmer  Delivery  Model  used  in  [Cunningham  (1993)]  for  which  a  set 
of  hydrodynamic  coefficients  and  geometric  properties  is  available.  Although 
the  main  results  and  conclusions  of  this  work  are  derived  for  a  submerged 
vehicle,  similar  techniques  can  be  applied  to  surface  ships  as  well. 
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II.  EQUATIONS  OF  MOTION 


A.  COORDINATE  SYSTEM 

In  our  analysis  we  are  going  to  use  a  moving  coordinate  system  {x,y,z), 
attached  on  the  vehicle.  The  origin  of  this  system  coincides  with  the  center 
of  buoyancy,  B.  The  x-axis  is  attached  to  the  longitudinal  plane  of  symmetry 
for  the  vehicle,  the  y-axis  is  positive  starboard,  and  the  z-axis  is  positive 
downwards.  All  main  symbols  used  in  the  development  of  the  equations  of 
motion  in  this  work  are  summarized  in  Table  1. 


B.  GENERAL  FORM  OF  THE  EQUATIONS  OF  MOTION 

The  equations  of  motion  for  a  submersible  vehicle  in  the  horizontal  plane 
are  written  as  follows: 

Sway  equation: 


m[v  +  Ur  -wp  +  xaipq  +  r)  -  yaip^  +  r^)  +  zaiqr  -  p)]  = 
Ypp  +  Yi-r  +  Ypqpq  +  Yqrqr  +  v  +  Yy^vw  + 

Ys^U^Sr  +  YyV  +  YpUp  H-  Yrll r  +  Yygvq  +  Y^p-wp  +  Y^rWr  + 


{W  —  B)  cos  OsiTKf)  — 

f  [CDyh{x){v  +  xrf  +  C'j5^?)(x)(u)  -  xq)^]  dx  . 

Ucf[X) 


(1) 
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Yaw  equation: 


Izzi"  +  {lyy  Ixx)PQ.  Ixy{p^  ~  9^)  ~  Iyz{pT  +  9)  +  Ixz{<l'<^  ~  P)  + 

m[xG{v  +  Ur  —  wp)  —  yG{U  —  vr  +  wq)]  = 

Npp  +  Nj-r  +  Npgpq  +  Ngrqr  +  N.ijUv  +  NyyjVw  + 

Ns^U'^Sr  +  Ni,v  +  NpUp  +  NrlJ  r  +  N^gvq  +  Nyjpwp  +  N^rwr  + 

{xqW  -  xbB)  cos  6*  siiK^  +  {yoW  -  ysB)  sin  9  +  U'^Np^op  - 
[  [<^Dj,^(a:)(v  +  a:r)^  +  C£.^6(a:)(w  -  X9)^l  .  (2) 

■’  Ucf{x) 

Roll  equation: 

IxxP  +  i^zz  ~  ^yy)Q'^  +  ^xy{pf'  ~  4)  ~  Iyz{<l^  ~  ~  Ixz{PQ  +  ^^ )  + 

mlyaiw  —  Uq  +  vp)  —  zg{v  +  Ur  —  wp)]  = 

Kpp  +  Kj-r  +  Kpgpq  +  Kg-^qr  +  KyU  v  +  Kyyjvw  + 

KpropU^  +  Kyi)  +  KpUp  +  KyU  r  +  KygVq  +  Kyjpwp  +  K  yjrWr  + 

(pqW  -  ysB)  cos  Bcoscp-  (zaW  -  zbB)  cos  9sin<f>  .  (3) 

The  rotational  velocity  equation  around  the  x-axis  is  simply, 

^  =  P ,  (4) 

while  Ucf  denotes  the  cross-flow  velocity, 

Ucf{x)  =  -\/{v  +  xr)2  +  {w  —  xq)^  . 
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b{x) 

local  beam  of  the  hull 

Cd 

quadratic  drag  coefficient 

Sr 

rudder  deflection 

h{x) 

local  height  of  hull 

lyy'i  Izz) 

vehicle  mass  moments  of  inertia  about  body  axes 

1 

cross  products  of  inertia 

vehicle  mass 

Euler  angles 

u 

constant  vehicle  speed  along  the  x-axis 

{u,v,w) 

translational  velocities  about  {x,y,z)  axes 

distances  along  the  three  body  axes 

(X,Y,Z) 

force  components  along  the  body  axes 

coordinates  of  the  center  of  gravity 

coordinates  of  the  center  of  buoyancy 

^nose 

fore  coordinate  of  vehicle  body 

^tail 

aft  coordinate  of  vehicle  body 

Table  1;  Nomenclature 


C.  SIMPLIFICATIONS 

Before  we  proceed  with  the  analysis,  we  must  simplify  the  above  equations 
of  motion  in  order  to  reflect  the  fact  that  we  are  analyzing  motion  in  the 
horizontal  plane  only.  The  simplifications  that  we  employ  are  the  following: 


•  Velocity,  ty,  and  acceleration,  ru,  in  the  ^-direction  are  zero. 

•  Acceleration  in  the  longitudinal  direction,  tt,  is  zero. 

•  Rotational  velocity,  and  acceleration,  g,  in  the  y— direction  are  zero. 

•  Pitch  angle,  6,  is  zero. 
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D.  SIMPLIFIED  EQUATIONS  OF  MOTION 


After  the  application  of  the  above  simplifying  assumptions,  the  equations 
of  motion  become: 

Sway  equation: 

(m  -  Yij)v  +  {mxQ  -  Yf)r  -  [mza  +  Yp)p  = 

YyUv  +  {YrU  -  mU)r  +  YpUp  +  ycp^  +  yar^  +  (W  -  B)  sin.^  +  Ys^U^Sr  - 

(5; 

Ucf{x) 

Yaw  equation: 

(mxG  -  Ni)v  +  -  Nf)r  -  {I^z  +  Np)p  = 

NvU  V  +  (NrU  —  mxoUy  +  NpUp  +  Ixyp"^  +  lyzpr  +  yovr  + 

{xgW  -  xbB)  sin  4>  +  Ns^U^Sr  +  U^Np^op  - 

r””Co.Mx)(»+x.)^%±^xdx^  (6) 

Ucf{x) 

Roll  equation: 

(  TTIZg  Kv^V  "t"  (/xz  +  (/xx  ~ 

K.uU  V  +  [KrU  +  mzGU)r  +  KpUp  —  IxyPf  —  lyz'r'^  —  Tnyc^p  + 

{^gW  -  ysB)  cos  (p  -  {zgW  -  zbB)  sin  p  +  U^Kp^op  ■  (7) 


Roll  rate: 


p  =  p  . 


(8) 


III.  LINEAR  ANALYSIS 


A.  LINEARIZATION 

The  simplified  equations  of  motion  can  be  written  in  matrix  form  as  follows; 

Ax  =  Bx  +  g(x),  (9) 

where  the  state  vector  x  is  defined  as, 

V 

r 
P 
<!> 

and  the  state  matrices  are, 

m  —  Yy  mxo  —  Yf 

TTIXq  N i)  ^ zz  ^ f 
TTiZQ  Ky  Ixz  K  r 
0  0 

and 

■  YyU  YrU  -  mU  YpU  0  ' 

^  NyU  -mxcU  +  NrU  NpU  0 

KyU  mzaU  +  KrU  KpU  0  ' 

.0  0  1  0  _ 

The  term  g(x)  contains  all  the  nonlinear  terms, 

91  =  +  yar^  +  {W  -  B)sm4)  +  Y^JU'^Sy  - 

r^nose 

Coy  h{x){v  +  xr)\v  +  xr\dx  , 

92  =  Ixyp"^  +  lyzpr  +  VGvr  +  {xqW  -  xbB)  sin  4>  +  Ns^U'^Sy 

rXnose 

^Dy  h{x){v  +  xr)\v  +  xr\x  dx  , 

9s  =  -Ixypr  -  lyzr'^ -myGvp  +  U'^Kprop  + 


(10) 

(11) 


—mzG  —  Yp  O' 

-Ixz  -Np  0 
Ixx  —  Kp  0  ’ 

0  1 
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(^yaW  —  vbB)  cos  (j)  —  {zgW  —  zbB)  sin  <f) 


(12) 

(13) 


=  0. 

We  want  to  linearize  the  nonlinear  terms  about  a  nominal  point 

XQ  =  bo,ro,Po,<?^o]^  =  0  . 

After  applying  Taylor  series  expansion  of  the  non-linear  terms  about  the  nom¬ 
inal  point  xq  and  keeping  only  the  linear  components,  the  linearized  equations 
of  motion  are  written  in  matrix  form  as: 

A'x  =  B'x  ,  (14) 

where 

A'  =  A  , 


and 

■  YvU 

YrU  —  mil 

YpU 

W  -  B 

B'  = 

N^U 

—mxaU  +  NrU 

NpU 

xqW  —  xbB 

K^U 

mzcU  +  KrU 

KpU 

—zqW  -t-  ZbB 

/  ^  A  \  m 

0 

0 

1 

0 

Equation  (14)  is  our  linearized  system.  Eigenvalue  analysis  for  this  system  is 
performed  in  the  next  section  in  order  to  assess  the  dynamic  stability  of  the 
vehicle. 


B.  LOSS  OF  STABILITY 

Stability  of  the  linearized  system  depends  on  the  location  of  the  four  eigen¬ 
values  of  the  system,  in  other  words  the  roots  of  the  characteristic  equation: 

det(B'  -  XA')  =  0  ,  (15) 
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or,  in  polynomial  form, 


-<4A  "h  B +  CA^  +  DX  E  =  0  .  (16) 

The  coefficients  of  the  polynomial  equation  (16)  are  given  by: 

^  =  (Ixx  —  Kp)[{m  —  Yy){Iz2  —  Nf)  —  {mxc  —  Yr){mxo  —  Ny)] 

+  {mzG  +  K.o)[{—Izz  +  Nr){mzG  +  Y^)  +  {Nj,  +  Ixz){mxG  -  Tr)] 

-  {Kr)[{-mxG  +  Ni){mzG  +  Fp)  +  {Np  +  7^^)(m  -  y^,)]  ,  (17) 

B  =  KvU[{-Izz  +  Nf){mzG +  Yp)  +  {Np  + Ixz){mxG -Yf)] 

—  {mzG  +  Ky)[{Izz  —  Nr){YpU)  —  {U  Nr  —  UmxG){'niZG  +  Yp) 

+  {Np  +  hz){UYr  -  Urn)  -  {NpU){mxG  -  F)] 

-  {hx  -  Kp)[{Tn  -  Yy){U Nr  -  UmxG)  +  {YyU){Izz  -  Nf) 

—  {mxG  —  Yr){NyU)  —  {mxG  —  Ny){UYr  -  Um)] 

—  KpU[{m  —  Yy){Izz  —  Nr)  -  {mxG  —  Yr)(mxG  -  Ny)] 

+  {Kr  -  hz)[{mxG  -  Ny)(YpU)  -  {NyU){mzG  +  Fp) 

+  (Np  +  hz){YyU)  -  {NpU){m  -  F^,)] 

—  (mzGU  +  KrU)[{—mxG  +  Ni,){mzG +  Yp) 

F  (-A^p  F  Ixz){'i^  Yi])\  (1^) 

C  =  [Ixx- Kp)[(YyU){UNr-mUxG)- {UYr-Um){NyU)] 

+  KpU[{m  -  Yy){UNr  -  mUxG)  F  (YyU){Izz  -  Nr) 

-  (mrrG  -  Yr){NyU)  -  {mxG  -  Ny){UYr  -  Um)] 

F  (mzGU  +  KrU)[{mxG  —  Ny){YpU)  —  {NyU){mzG  +  Fp) 

F  (iVp  +  7,,)(F„17)-(iVpt/)(m-Fi,)] 
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+  {zgW  -  ZBB)[{m  -  Yy){lz2:  -  Nr)  -  {mXG  -  Yr){mXG  -  Nv)] 

-  {mzG  +  Ki,)[{xBB  -  XGW){mxG  -  Yr)  -  {U Nr  -  UmxG){YpU) 

+  {UYr  -  Um)(NpU)] 

—  KyU[{Izz  —  Nr){YpU)  —  {U  Nr  —  UmxG)(mzG  +  Yp) 

+  {Np  +  Ixz){UYr  -  Um)  -  {NpU){mxG  -  Yr)] 

+  {Kr  —  Ixz)[{xbB  -  XGW){m  -  Yy) 

-  {N,U){YpU)  +  iNpU){YM)]  ,  (19) 

D  —  {mzGU  +  KrU)[{xBB  —  XGW){m —  Yp) 

-  {N,U){YpU)  +  {NpU){Y,U)] 

-  {Kr  -  1xz){xbB  -  XGW){YyU) 

+  {mzG  +  Kp){xBB -XGW){UYr-Um) 

—  KyU[{xBB  —  XGW){mXG  —  Yr) 

-  {UNr-mUxG){YpU)  +  {UYr-Um){NpU)] 

-  KpU[{Y^U){UNr  -  mUxG)  -  {UYr  -  Um){N^U)] 

-  {zG^  -  ZBB)\{m  -  Yy){UNr  -  mUxG)  +  (YvU){Izz  -  Nr) 

-  {mxG  -  Yr){N^U)  -  [mxG  -  Np){UYr  -  Um)]  ,  (20) 

E  =  {zgW  -  ZBB)[{Y^U){UNr-mUxG)- {UYr-Um){N^U)] 

+  {K^U){xbB  -  XGW){UYr  -  Um) 

-  {mzGU  +  KrU){xBB  -  xgW){Y^U)  (21) 

We  can  examine  the  stability  of  the  system  by  utilizing  Routh’s  criterion. 
Application  of  this  criterion  to  the  characteristic  equation  (16),  reveals  that 
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the  following  two  conditions  must  be  satisfied  in  order  to  ensure  that  all  roots 
of  (16)  have  negative  real  parts: 


BCD  -  AD^  -  EB^  >  0  , 

(22) 

E  >0. 

(23) 

If  E  is  less  than  zero,  one  real  root  of  (16)  becomes  positive  and  the  sytem  will 
become  unstable  in  a  divergent  manner  (Guckenheimer  and  Holmes,  1983). 
This  is  the  case  of  a  directionally  unstable  ship  which  is  well  known  in  the 
literature  (Clayton  and  Bishop,  1982).  If,  however,  condition  (22)  is  violated, 
the  system  will  exhibit  an  oscillatory  motion  due  to  the  presence  of  complex 
conjugate  roots  with  positive  real  parts.  This  form  of  instability  is  caused  by 
the  coupling  of  roll  into  sway  and  yaw  and  is  further  analyzed  in  this  work. 

In  order  to  compute  the  limiting  case  of  loss  of  stability,  we  consider  equa¬ 
tion  (24), 

BCD  -  -  EB^  =  0  .  (24) 

The  result  of  this  equation  will  produce  the  limiting  value  of  zq  as  a  function 
of  XG  for  loss  of  stability.  A  curve  of  this  functional  form. 


ZG  =  f{xG)  , 

will  be  our  locus  of  loss  of  stability.  After  some  algebra,  we  can  express  the 
coeflficients  of  equation  (24)  in  the  following  form: 

A  =  Ai4  +  A2ZG  +  As  ,  (25) 
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where 


Ai  =  -  Nr) 

A2  =  —^¥^{122  —  Nr)  —  mKij{Izz  -  Nr)  +  m{Np  +  Ixz){mxo  —  Yr) 

+  m{Kf  -  Ixz){mxG  -  Ni;) 

M  =  (hx  -  Kp){m -Yy){lzz- Nr)  -  {Ixx- Kp){mxG -Yr){mxG  -  Ny) 

-  KyYp{Ixx  -  Nr)  +  Kv{Np  +  Ixz){'mXG  -  Yr) 

+  Yp{Kf  -  Ixz){mxG  -  Ni)  -  (Kr  -  Ixz)iNp  +  Ixz)im  -  Yi) 

B  =  BiZq  +  B2ZG  +  ,  (26) 

where 

Bi  =  m^{UNr  —  Utuxg)  +  m‘^U(mxG  —  Ni) 

B2  =  -m{KyU){Ixz- Nr)  -  m{Ixz- Nr){YpU) 

+  mYp{UNr  —  UmxG)  +  TnKi^UNr  —  Uttixg) 

—  m{Np  +  Ixz){UYr  —  Um)  +  m{NpU){mxG  —  ^r) 

—  m{Kr  —  Ixz){NvU)  +  mUYp{mxG  —  Ni) 

—  mU{Np  + Ixz)i‘m —  Yi) +  mUKr(mxG  —  Ni) 

Bz  =  -Yp{KMhz-Nr)  +  {KMNi,  +  Ixz){mxG-Yr) 

-  Ki{Ixz-Nr){YpU)AKiY0Nr-UmxG) 

—  Ki{Np  +  Ixz){UYr  —  Um)  +  Ki{NpU){mxG  —  Yj) 

-  {I XX  -  Kp){m  -  Yi){UNr  -  UmxG)  -  {Ixx  -  Kp){YyU){Ixz  ~  Nr) 

+  {Ixx  -  Kp){mxG  -  Yr){NyU)  +  {Ixx  -  Kp){mxG  -  Ni){UYr  -  Um) 
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-  iKpU)im  -  n)(/,,  -  Nr)  +  (KpU)imxG  -  Yr}{mxG  -  N^) 

+  {Kr  -  Ixz){mxG  -  Ny){YpU)  -  Yp{Kr  -  Ixz){NyU) 

+  {Np  +  Ixz){Y^U){Kr  -  Ixz)  -  {Kr  -  Ixz){NpU){m  -  Yi) 

+  UKrYp{mxG  —  Ny) 

C  =  Ciz%  +  C2ZG  +  Cz  ,  (27) 

where 

Cl  =  -m?U{NyU) 

C2  =  mU {mxG  -  Ny){YpU)  -  mUYp{NyU)  -  mUKr{NyU) 

+  mU{Np  + Ixz)iYyU) —  mU{NpU){m —  Yi,) 

+  W{m  -  Yy){Izz  -  Nr)  - W {mxG  -  Yr){mxG  -  Ny) 

-  m{XBB  —  XGW){mxG  -  Yf)  +  m{UNr  —  UmxG){YpU) 

-  m{UYr  -  Um){NpU)  +  m{KyU){UNr  -  Utuxg) 

Cz  =  {Ixx  -  Kp){YyU){UNr  -  UmxG)  -  (hx  -  Kp){UYr  -  Um){NyU) 

+  {KpU){m  -  Yy){UNr  -  Uttixg)  +  {KpU){YyU){Izz  -  Nj) 

+  {KpU){mxG  -  Yr){NyU)  -  {KpU){mxG  -  Ny){UYr  -  Um) 

+  UKr{mxG  -  Ny){YpU)  -  UKrYp{NyU) 

+  UKr{Np  +  Ixz){YyU)  -  UKr{NpU){m  -  Yy) 

-  Ky{XBB  -  XGW){mxG  -  Yr)  +  Ky{UNr  -  UmxG){YpU) 

-  Ky{UYr  -  Um){NpU)  -  {KyU){Izz  -  Nr){YpU) 

+  Yp{KyU){U  Nr  —  U  ttlxg)  —  {KyU )  {Np  +  Ixz)  {UYr  —  U  m) 

+  {KyU){NpU){mxG  -  Yr)  +  {Kr  -  Ixz){XbB  -  XGW){m  -  Yy) 
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(28) 


-  (Kf  -  h.){KU){Y,U)  +  (K,  -  h.)(N,U){Y„U) 

D  =  Dizg  +  D2  , 

where 

D\  =  mU{xBB  -  XGW){'m —  Yi,)  -  mU{NyU){YpU) 

+  mU{NpU){Y^U)+m{xBB  XGW){UYr  -  Um) 

-  W{m-  n)(C/iV,  -  UmxG)  -  W{Y^U){I,^  -  Nr) 

+  W (mxG  -  Yf){NvU)  +  W {mxG  -  Ny){UYr  -  Um) 

D2  =  UKr{xBB-XGW){m-Yi^-UKriNyU){YpU) 

+  U  Kr(NpU){YvU)  —  (Kr  —  Ixz)ixBB  —  XgW)(Y.uU) 

+  Ky{xBB  —  XGW){UYr  —  Um)  —  {KyU){xBB  —  XGW){mxG  —  Yj) 

+  {KyU){UNr  -  UmxG){YpU)  -  {KyU){UYr  -  Um){NpU) 

-  {KpU){YyU){UNr  -  UmxG)  +  {KpU){UKr  -  Um){NyU) 

and 

E  =  EiZG  +  E2  ,  (29) 

where 

El  =  W{YyU){UNr-UmxG)-W{UYr-Um){NyU) 

-  mU{xBB -XGW){YyU) 

E2  =  {KyU){xBB -XGW){UYr-Um)-UKr{xBB -XGW){YyU) 

If  we  apply  the  stability  criterion  (24)  utilizing  expressions  (25)  through 
(29),  we  get  a  fifth  order  polynomial  equation  in  the  metacentric  height  zg  of 
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the  following  form, 


+  F^zq  +  Fzz%  +  F2ZQ  +  Fizq  +  Fo  =  0  ,  (30) 

where, 

Fo  =  B3CZD2  -  AsDl  -  E2BI  , 

Fi  =  FsCsDi  +  (F3C2  +  F2C3)D2  -  2F2B2B3  -  - 

2D2D,A2  -  dIA2  , 

F2  =  -F2(F|  +  2F1F3)  -  2F1F253  +  (FsCa  +  F2C3)r>i  + 

(53^1  +  B2C2  +  BiC'3)r>2  -  F^^3  -  2D2D1A2  -  D^Ai  , 

F3  =  -DfA2-2D2DiAi-2E2BiB2- Ei{B^  +  2BiBz)  + 

DiiB^Ci  +  B2C2  +  BiCz)  +  D2{B2Ci  +  B1C2)  , 

F4  =  Fi(B2C'i  +  B1C2)  +  B1C1D2  -  D^Ai  -  E2BI  -  2E1B1B2  , 

Fg  =  -EiBf  +  BiCiDi  . 

Using  the  corresponding  Matlab  program  shown  in  Appendix  A,  we  can  solve 
equation  (30)  using  a  typical  constant  forward  velocity  U  —  5  ft/sec.  The  value 
of  xa  is  given  in  the  range  from  0  to  2  ft.  In  this  way  we  get  five  solutions  in 
ZQ  for  a  given  value  of  xq.  Investigating  these  solutions,  we  can  see  that  only 
one  satisfies  the  second  criterion  (23).  This  solution  corresponds  to  the  locus 
of  loss  of  stability.  For  values  of  zq  greater  than  the  critical  value,  the  system 
is  stable,  and  for  values  less  than  its  critical  value,  the  system  is  unstable,  see 
Figure  1.  The  calculations  are  repeated  for  a  range  of  forward  speeds  U  from 
2  to  8  ft/sec.  The  results  are  shown  parametrically  in  Figure  2.  We  observe 
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an  upwards  movement  of  the  curves  as  the  speed  is  increased.  This  means 
that  the  system  experiences  a  tendency  to  become  less  stable  at  higher  speeds 
which,  in  turn,  calls  for  higher  metacentric  heights  to  ensure  stability. 
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IV.  NONLINEAR  ANALYSIS 


A.  INTRODUCTION 

From  the  linearized  analysis  on  the  loss  of  stability  that  was  done  in  the 
previous  chapter,  we  can  see  that  as  a  specific  parameter  of  the  system,  such  as 
(xQ,  zg),  is  varied,  it  is  possible  to  pass  from  a  region  of  stability  to  a  region  of 
instability.  This  case  of  loss  of  stability  is  associated  with  one  pair  of  complex 
conjugate  eigenvalues  of  the  system  crossing  the  imaginary  axis.  This  loss  of 
stability  is  usually  accompanied  by  self  sustained  oscillations,  and  it  is  called 
Hopf  Bifurcation. 

There  are  two  cases  of  Hopf  Bifurcations,  supercritical  and  subcritical.  In 
the  supercritical  case,  limit  cycles  are  created  just  after  the  loss  of  stability. 
A  limit  cycle  is  a  constant  amplitude  oscillatory  motion  of  our  system.  This 
amplitude  is  usually  larger  as  we  move  away  from  the  bifurcation  point.  When 
the  limit  cycles  have  small  amplitudes  the  situation  is  not  very  critical  for  our 
vehicle  and  is  not  much  different  from  stable  conditions.  In  the  subcritical 
case  we  may  have  convergence  to  limit  cycles,  even  before  our  system  looses 
its  stability.  Furthermore,  the  limit  cycle  amplitudes  are  considerably  higher. 

The  analysis  that  will  follow  will  be  performed  in  order  to  verify  the  exis¬ 
tence  of  the  limit  cycles  and  to  find  out  which  case  of  Hopf  Bifurcations  we 
have  in  our  model.  We  will  also  examine  the  stability  of  the  resulting  limit 
cycles.  This  analysis  is  necessary  because  we  can  predict  the  behavior  of  a 
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vehicle  when  some  of  its  parameters  change  and  we  have  operation  in  the 
proximity  of  a  bifurcation  point.  The  results  of  the  non  linear  analysis  will 
be  verified  by  a  numerical  simulation  of  the  system’s  motion,  both  above  and 
below  the  bifurcation  point. 


B.  THIRD  ORDER  EXPANSIONS 

From  the  linearization  procedure  we  see  that  our  system  is  written  in  the 
form  of  matrix  equation  (14),  where  we  have  ignored  the  non  linear  terms. 
If  we  take  into  account  the  non  linear  terms  up  to  third  order,  equation  (14) 
becomes: 

Ax  =  B'x  +  g(x)  ,  (31) 

where, 

9z{x) 

9a{x)  _ 

Keeping  terms  up  to  third  order  the  vector  of  non  linear  terms  can  be  written: 

g{x)  =  +  g^^\x)  +  c(x)  ,  (32) 

where  g^‘^\x)  contains  the  second  order  non  linear  terms,  g^^\x)  contains  the 
third  order  non  linear  terms,  and  c(x)  contains  the  constant  terms.  The  cross 
flow  integrals  can  be  written  as  follows: 

rXnose 

I  h{x){y  xr)\v  +  xr\dx 

r^nose 

Ir  =  Coy  /  h{x){v  +  xr)\v  +  xr\x  dx 


g(x)  = 
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(33) 


The  second  order  non  linear  terms  are: 

9?\x) 
92\x) 
9Pix) 

94\^) 

where 


9i  =  vgP^  +  yar"^  - 

92^  =  hyV^  +  lyzpr  +  yQVr  - 

/Q\  _ 

53  =  -hypr  -  lyzT  -  myavr  -  {yaW  -  yBB)(j)‘^ 

=  0. 


The  third  order  non  linear  terms  are: 


where, 


g^^^(x)  = 


5^(3:) 
92  \x) 
9z\x) 

94\x) 


5^  = 

0 

92^  =  -4^^  -  ^{xgW  -  xbB)4>^ 
9^  =  ^izcW  -  zbB)<I>^ 

5f  =  0. 


The  constant  terms  are: 


c(x) 


■  ci{x)  ■ 
C2ix) 
czix) 
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where 


ci{x)  =  YsrU\ 

C2(x)  =  NsM\  +  U^Nprop 

C3{x)  =  U'^Kprop  +  VgW  -  VbB 

c^{x)  =  0  . 

In  order  to  get  the  second  and  third  order  non  linear  terms  of  the  cross  flow 
integrals,  we  must  expand  in  Taylor  series  a  function  of  the  form: 


fio  =  m 


(34) 


about  a  nominal  point  ^o- 


Jill  =  folJol  +  2|{ol(f  -  Jo)  +  sign({o)(J  -  Jo)“  + 


(35) 


The  sign  function  in  equation  (33)  can  be  approximated  by: 

sign(^o)  =  linitanh(— )  , 

7-*0  'y 

where  the  approximation  gets  better  as  7  gets  smaller. 

If  we  choose  ^0  =  0  as  our  nominal  point,  equation  (35)  becomes: 


(36) 


JIJI  =  J-J® 

67 


In  our  case  ^  is  (v  +  xr)  so  we  have: 


(37) 


or 


(v  +  xr)\v  +  xr\  =  — (v  +  xr)^ 
67 


{v  +  xr)\v  +  xr\  —  ;;^(v^  +  +  Zv^xr  +  3xVt>) 

07 


(38) 


(39) 
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Using  equation  (39)  the  cross  flow  integrals  become: 

{Eqv^  +  SEiv'^r  +  3E2vr^  +  E^r^)  (40) 

It  =  -^iEiv^  +  3E2v\  +  3E3vr^  +  E4r^)  (41) 

where, 

r^nose 

Ei=  x^h{x)dx  ,i  -  0,1, 2, 3, 4  (42) 

By  using  this  approximation  we  see  that  the  expansion  of  the  cross  flow  inte¬ 
grals  give  only  third  order  terms,  so  we  have: 


/f  =  =  0 


(43) 


Now  we  want  to  write  our  matrix  equation  in  state  space  form,  so  we  have  to 

find  the  inverse  of  the  system  matrix: 

£11  £12.  £12.  n  1 

D  D  D  ^ 

£21  022  023  n 

D  D  D  ^ 

Ml  0.32  0.33  A  ’ 

D  D  D  ^ 

0  0  0  1 


where, 


{Izz  IIr){Ixx  Kp)  {Ixz  ~  Kf){  —  Ixz  ~  Up) 

«12  =  (Xr  —  mXG){lxx  —  Kp)  -|-  (Ixz  —  Kr)(—mZG  —  Yp) 
ai3  =  (mxG  -  Yr){-Ixz  -  Np)  -  {I^^z  -  Nr){-mzG  -  Yp) 

021  =  -{mxG- Ny){I^^- Kp)  +  {-mzG- Ki,){-l^x- Np) 

022  =  {m-  Yi){lxx  -  Kp)  -  (-mzG  -  Ky){-mzG  -  Yp) 

023  =  -Yy){-Ixz  -  Np)  +  (mxG  -  Ni,)(-mzG -Yp) 
031  =  (mxG  -  Ny){Ixz  -  Kr)  -  {-mzG  -  Ky){I;zz  -  Nf) 
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032 


—  —  (m  —  Y.u){Ixz  —  Kr)  +  {-mzG  —  Ky){mxG  —  Yr) 
033  =  {m -Yij){Izz- Nf)  -  {mxG  -  Ny){mxG-Yr) 

and 

D  =  {m-Y^{hx-Nr){Ixx-Kp) 

-  {m-Yi,){Ixz- Kr){-Ixz- Np) 

{tYIXg  ^v'}(jXLXg  Yx^(^Ixx 
+  {mxG  -  Ni)(Ixz  -  Kf)i-mzG  -  Yp) 

+  {-mzG  -  Kv){mxG -Yr){-Ixz  -  Np) 

-  {-mzG  -  Kij){Ixz  -  Nr){-mzG  -  Yp) 

If  we  multiply  equation  (31)  by  {A'y^  from  the  left  side  we  get: 


{A')-^A'i  =  {A'y^B'x  +  {A')~'^g{x) 


or 


X  =  Fx  +  G{x) 


where, 


F  =  {A')-^B\ 
G(x)  =  (A')-^g(x) 


Then  matrix  F  is  a  4  x  4  matrix  defined  as  follows; 


r  £ll  £12  £12,  £k  1 


D  D  D  D 


0  0  10 


(44) 


(45) 


(46) 

(47) 
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where 


Fn  =  an{Y^U)  +  aniKU)  +  an{KyU) 

Fi2  —  O’liiXr  —  fn)U  +  ai2{Nr  —  mxG)U  +  aiz{Kr  +  mzG)U 
Fiz  =  auiXpU)  +  auiNpU)  +  aiz{KpU) 

Fi4  =  axxiW  —  B)  +  aizixoW  —  xbB)  +  aiz{—ZQW  +  zbB) 

F21  —  0'2\{YvU^  +  a22{^N.JJ^ a2z{^K.yU) 

F22  =  o,2i{Yr  —  m)U  +  a22{Nr  —  TnxG)U  +  a23{Kr  +  mzG)U 

F23  =  o.2i{YpU)  +  a22{NpU)  +  a23{KpU) 

F24  =  (^2iiW' —  B)  +  a22{xGW  —  xbB)  +  a23i—ZGW  +  zbB) 

F31  =  asiiY^U)  +  a32{NyU)  +  assiK^U) 

F32  =  0,3l{yr  ~ '^)u  +  az2{Nr  — 'mXG)U  +  a33{Kr  +  mZG)U 
F33  =  o,zi{YpU)  +  a32{NpU)  +  az3{KpU) 

-P’34  =  <^3iiW  ~  B) -\- az2{xGW  —  xbB)  +  az3{—ZGW  +  zbB) 

From  equations  (45), (47)  we  see  that: 


where, 


G{x)  =  {A')-^g{x)  = 


■  G'i(x)  ■ 
G2{x) 
G3{x) 

.  G4{x)  . 


Gi(a:)  =  -^^i(x)  + -^^2(2:)  +  ^53(2:) 

G2{.x)  =  -^g-^[x) +  ^g2{x) +  ^g2,{x) 

G3{x)  =  ^g^[x)  +  ^g2{x)  +  ^g^{x) 

Gi{x)  =  0 
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Equation  (45)  can  also  be  written  in  the  following  way: 

i  =  Fi+GW(i)  +  G®(x)  +  i: 


(48) 


where, 


G^^\x)  = 


G^^\x)  = 


and 


k  = 


Gf(x) 

Gf\x) 

G?\x) 

Gf(x) 

G?(x) 

G^Hx) 

Gf\x) 

Gf\x) 


ki 

k2 

h 

ki 


Each  element  of  the  above  non  linear  terms  vectors  can  be  written  as  follows: 


G?\x) 

GfHx) 

afHx) 

Gf(x) 


“11  (2).  X  ,  “12  (2)/  N  ,  “13  (2)/  V 

-^9i  '{x)  +  —9^  \x)  +  — \x) 

“21  (2)/  X  ,  “22  (2)/  X  ,  “23  (2)/  x 

-^9i  \x)  +  —9^2  +  -^93 

“31  (2)/  X  ,  “32  (2),  X  ,  “33  (2),  x 

-^9l  \x)  +  —9^2  (^)  +  -^^3  (^) 

0 


and 


Gfix) 


!?i 
Gf\x) 
afHx) 
Gf(x) 


“11  (3)/  X  ,  “12  (3)/  X  ,  “13  (3)/  X 

~^9i  {^)  +  ^^92  (x) 

“21  (3)/  X  ,  “22  (3)/  X  ,  “23  (3)/  x 

■^91  \x)  +  —9^  >{x)  +  —9l  ^(x) 

“31  (3)/  X  ,  “32  (3)/  X  ,  “33  (3)/  x 

-^9l  +  -^92  +  —9\  >{x) 


0 
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Finally,  the  constant  terms  are: 


On 

+ 

Oi2 

Ol3 

D 

Cl 

D 

C2 

+ 

D 

021 

+ 

022 

+ 

023 

D 

Cl 

D 

C2 

D 

^31 

+ 

032 

+ 

033 

D 

Cl 

D 

C2 

D  """ 

A:4  =  0. 


C.  COORDINATE  TRANSFORMATIONS 

From  equations  (45)  and  (48)  it  is  obvious  that  the  stability  of  our  system 
depends  on  the  eigenvalues  of  matrix  F.  Since  we  want  to  investigate  the 
behavior  of  our  system  oround  the  Hopf  Bifurcation  point,  it  is  useful  to  bring 
our  system  in  its  normal  coordinate  form.  This  can  be  done  by  applying  a 
transformation  of  the  coordinate  system,  using  as  transformation  matrix  T, 
the  modal  matrix  of  eigenvectors  of  F,  evaluated  at  a  critical  point.  This 
critical  point  will  be  a  pair  of  zg  and  xq  values,  that  belong  to  the  critical  line 
of  loss  of  stability.  The  applied  transformation  will  be  as  follows: 


X  =  Tz 


(49) 


or. 


2  =  T-^x 

Then  equation  (48)  can  be  written  as  follows; 


(50) 


Tz  =  FTz  +  {Tz)  +  {Tz)  Fk  (51) 
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or, 

i  =  T~^FTz  +  T-'^G^‘^\Tz)  +  T~^G^^\Tz)  +  T~'^k  (52) 

In  this  new  coordinate  system,  the  system’s  matrix  is  T-'^FT,  and  at  the  Hopf 
Bifurcation  point  can  be  written  as  follows: 

0  —ojQ  0  0 

rp— iTTinn  ^0  0  0 

^  0  0  P,  0  ' 

.0  0  0  P2  . 

Where  ojq  is  the  imaginary  part  of  the  critical  pair  of  eigenvalues,  and  the 
remaining  eigenvalues  Pi  ,  P2  are  negative.  For  values  close  to  the  Hopf  Bi¬ 
furcation  line,  we  can  write  the  system  matrix  as  follows: 

off.  —{uiQ  +  uj'e^  0  0 

rji — ^FT (^0  -f-  aJ  e)  a  €  0  0 

“0  0  (Pi  +  Pi'e)  0 

0  0  0  (P2  +  P2^)  . 

where:  e  =  Difference  from  the  critical  point  {zq  —  zq^). 

—  Derivative  of  the  real  part  of  the  critical  eigenvalue  with  respect  to  e. 

=  Derivative  of  the  imaginary  part  of  the  critical  eigenvaluewith  respect 

to  €. 

P[  =  Derivative  of  Pi  with  respect  to  e. 

P2  =  Derivative  of  P2  with  respect  to  e. 

D.  CENTER  MANIFOLD  EXPANSIONS 

By  using  the  coordinate  transformation  described  in  the  previous  chapter. 
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and  from  equation  (50),  we  see  that  we  have  a  new  variables  vector,  which  is: 

"  zi  ' 

^3 

.  ^4  . 

The  first  two  coordinates  zi  ,  Z2,  are  the  critical  coordinates  and  correspond 
to  a  pair  of  complex  conjugate  eigenvalues.  The  remaining  two  coordinates 
^3,  Z4,  are  the  stable  coordinates  and  they  always  correspond  to  eigenvalues 
that  are  negative.  The  center  manifold  theory  predicts  that  the  relationship 
between  the  critical  coordinates  zi,  Z2,  and  the  stable  coordinates  z^,  Z4,  is  at 
least  of  quadratic  order.  After  this  assumption,  the  two  stable  coordinates  can 
be  written  as  follows: 


2  2 

2:3  =  h\\Z^  +  hi2ZiZ2  +  h22^2 

(53) 

2  0 

Z4  =  SiiZ-^  +  S12Z1Z2  +  S22Z2 

(54) 

The  coefficients  hij,  Sij,  of  equations  (53),  (54),  need  to  be  determined.  By 
differentiating  equations  (53),  (54),  we  get  the  following: 

is  =  2hiizizi  +  hi2{ziZ2  +  Z1Z2)  +  2h22Z2Z2  (55) 

Z4  =  2suzizi  +  si2{ziZ2  +  Z1Z2)  +  2S22Z2Z2  (56) 

From  equation  (52)  and  if  we  ignore  the  higher  order  terms,  we  take: 

=  —u)oZ2  (57) 

22  =  UQZi  (58) 

If  we  substitute  equations  (57),  (58),  into  equations  (55),  (56),  we  get  the 
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following: 


^3  —  hi2U>o^i  +  2(/i22  -  hn)‘^o^iZ2  —  hi20iJoZ2  (59) 

Z4  =  512^0^1  +  2(s22  -  Sn)uioZiZ2  -  Si2i^o^2  (60) 

The  third  and  the  fourth  equations  of  matrix  equation  (52)  can  be  written  in 
the  following  way: 


4  =  ^^i^3  +  [T-'G'(2)(r2:)](3,3)+r-'A:(3.3)  (61) 

i4  =  P2Z4  +  [T-^G^^\Tz)]^4A)  +  T-^h^,4)  (62) 

In  equations  (61),  (62),  we  kept  terms  up  to  second  order. 

The  transformation  matrix  T  and  it’s  inverse  T~'^  are  4x4  matrices,  and 
their  elements  can  be  denoted  by: 


T  =  [mij],  T-^  =  [m^l  i,  j  =  1,  2, 3, 4 


From  equation  (52)  we  have: 


where, 


r-i(?(2)(rz) 


di 

d2 

dz 

d4 


di  =  niiGf\Tz)  +  ni.,a'i\Tz)+nuGf\Tz) 
dl  =  n2iGf(Tz)  +  n^G''A(Tz)+n^Gf\Tz) 
ds  =  n3iGA(Tz)  +  n3^G^A(Tz)  +  n33Gf\Tz) 
d,  =  n„G^A{Tz)  +  n„af\Tz)  +  nt3Gf'’{Tz) 


(63) 


(64) 

(65) 

(66) 

(67) 
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Also  from  coordinate  transformation,  the  relationship  between  the  old  and 
new  coordinates  is  as  follows; 


V  =  miizi  +  mi2Z2  +  rriizzs  +  muZ4  (68) 

r  =  m2izi  +  m22Z2  +  ^12323  +  m24Z4  (69) 

p  —  msizi  +  mz2Z2  +  m33Z3  +  mz^z^  (70) 

4>  =  m4i2i  +  7714222  +  TTl^zZz  +  7774424  (71) 


Finally  if  we  substitute  equations  (68),  (69),  (70),  (71),  and  the  expressions  for 
G\,  G2,  Gz,  G4  into  equations  (64),  (65),  (66),  (67),  we  get  the  final  expressions 
for  the  coefficients  di\ 


di 

= 

+  h6^1^2  +  ^17^2) 

+ 

'^12(125^1  +  l26^lZ2  +  I27Z2) 

+ 

^13(^35^1  h6^lZ2  +  ^7^2) 

(72) 

d2 

^21(^15^1  +  I17Z2) 

+ 

^22(^25^1  +  I2QZ1Z2  +  I27Z2) 

+ 

'^2z{hbz\  +  heZiZ2  +  I37Z2) 

(73) 

dz 

= 

+  Wq^IZ2  +  I17Z2) 

+ 

^32(^25^1  +  ^26^1'2^2  +  I27Z2) 

+ 

^33(^35^1  +  h^ziZ2  +  137  zl) 

(74) 

^4 

= 

'^4l{h5Zi  +  ll%ZiZ2  +  I17Z2) 

+ 

^42(^25^1  +  heZiZ2  +  ^27^2^2) 

+ 

^43(^35^2^1  +  h6ZiZ2  +  I37Z2) 

(75) 
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Coefficients  lij,i  =  1,2,3  j  =  5,6,7  in  equations  (72),  (73),  (74),  (75),  are  as 
follows: 

^1,5  = 

+  4-  +  yGm2imii) 

-  —  (7xym3im2i  +  Iyzm2i  +  myomumsi 

+  yoWmli-yBBrnl^)  (76) 

^11 

h,6  =  -^yG(2m3im32  +  2m2i  +  77222) 

+  ^[24^7713177732  +  Iyz{mzim22  +  77732?7l2l) 

+  yG(^2l777l2  +  ?7722777ii)] 

-  (7773177122  +  7773277721)  +  2/yz7772  l77722 

+  777yG(777ii77732  +  777i27773i)  +  2{yGW  -  ^55)7774177742]  (77) 

h,7  =  ^yG(77732  +  77722) 

+  "^(4y77732  +  7yz7773277722  +  7/G7772277712) 

“  ^[4y 7773277722  +  4^777  22  +  777yc?777i277732  +  (yGl^  -  y55)77742]  (78) 

45  =  ^yG(7773i  +  77731) 

+  ^(4y  777|i  +  /yz7773i77721  +  yG7772l777ii ) 

®23/r  ,  r  2  , 

-  —(7xy 7773177721  +  /yz7772i  +  7777/(77771177731 

+  yGl^"^4l  -  yB5777^l)  (79) 

h,6  =  ^yG(27773i77732  +  277721  +  "722) 

^22 

+  -^[24y"i31^32  +  7y^  (7773177722  +  7773277721) 

+  yG("72l777l2  +  7772277711)] 
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-  -^[Ixy{'mz\m22  +  m32m2i)  +  2Iyzm2im22 

+  myG{miirnz2  +  'rrii2rn2i) +  2{yGW  -  yBB)m4im42]  (80) 

7  021  ,  2  ,  2  \ 

<2,7  =  -^yG{m^2  +  "^22) 

^22  /  r  2  T 

+  -^{^xy'm^2  +  Iyz'mz2m22  +  yG?7l2270i2) 

^23  r  9 

-  —[-^ij/Tn 3277122  +  Iyz'ni22  +  my 01^^12^22  +  (t/G^^  —  yS-B)m|2]  (81) 

,  O31  /  2  I  2  \ 

<3,5  =  -^ycym^i  +  mji) 

(232  .  2 

+  -^(lx3/"^31  +  4^"^3l”^21  +  2/0^2177711) 

^33  /  r  ,  r  2 

-  —(-<13/7773177721  +  -<yz7772i  +  777yG777ii7773i 

+  yell" Will  -  yB^TTl^l)  (82) 

^31 

<3,6  =  -^2/0(27773177732  +  277721  +  W722) 

+  -^[2-la:3/7773i77732  +  4z(7773i  77722  +  W732W72l) 

+  2/g(w7217T7i2  +  77722777ii)] 

-  -^[-?a;3/(wi3lWi22  +  Wi32W72l)  +  24^7772177722 

+  777yG(777ii77732  +  777i27773i)  +  2{yGW  -  yB5)7774i77742]  (83) 

h,i  =  ^2/g(w7|2  +  777^2) 

CI32  >  2 

+  'd  (4j,77732  +  4^7773277722  +  yGW722W7l2) 

-  -^[43/Wi32W722  +  42^^22  +  Wi2/GWll2Wi32  +  (2/g11^  -  2/B-B)777l2]  (84) 
Then  equations  (61),  (62),  can  be  written  as  follows: 

4  =  -<^1-23  +  C<3  +  63  (85) 

Z4  =  P2Z4  +  c<4  +  64  (80) 
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where  63, 64  come  from  the  constant  terms  as: 


and 


T-^k  = 


ei 

62 

63 

64 


ei  =  ^1^11  +  ^2^12  +  ^3^13 

62  =  kin2i  +  k2n22  +  k'in2z 

63  =  kiuzi  +  k2nz2  +  kzuzz 

64  —  kin^i  +  /c2n,42  +  kzti4z 

Using  equations  (53),  (54),  we  can  write  equations  (85),  (86),  as  follows: 

zz  =  Pi{hnZi  +  hi2ZiZ2  +  h22Z2)  +  dz  +  ez  (87) 

Z4  =  P2{suzf  +  S12Z1Z2  +  S22Z2)  +  ^4  +  64  (88) 

And  using  equations  (72),  (73),  (74),  (75),  they  become: 

h  =  (Pih  11  +  +  7^32^25  +  ^33^35)-^i 

+  {Pihi2  +  nzihz  +  ^32^26  +  '<Z2zh&)ziZ2 

+  {P\h22  +  nzilu  +  ^32^27  +  ^33^37)'Z2  +  63  (89) 

^4  =  (-^2611  +  ’^31^15  +  ^^32^25  +  ‘nzzhz)^! 

+  (-P2612  +  U3i^i6  +  nz2h6  +  nzzlzz)ziZ2 

+  {P2S22  +  ^31^17  +  ^32^27  +  ^33^37)22  +  64  (90) 

Comparing  the  coefficients  of  equations  (89),  (90),  with  the  coefficients  of 
equations  (59),  (60),  we  get: 

— Pihii  +  LO^hiz  =  nzili5  +  7732/25  +  ^^33/35 
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2ujQhii  —  P\h\2  +  2u)Qh22  =  nsi/ie  +  n32h6  + 


UJQhi2  P\h22  —  n^lln  +  ^32^27  +  ^33^37 

Solution  of  the  above  3x3  linear  system  of  equations,  gives  us  coefficients  /in, 
^12)  h22- 

Also  in  the  same  way: 

~-P2Sii  +  UI0S12  =  7i4i/i5  + 1142/25  + 1143/35 
— 2a>oSii  —  P2S12  +  2a;oS22  =  U41/16  +  1142/26  +  1143/36 
—010512  —  P2S22  =  Il4l/l7  +  1142/27  +  1143/37 

Solution  of  the  above  3x3  linear  system  of  equations,  gives  us  coefficients  sn, 
512  )  522- 

E.  AVERAGING 

In  this  part  of  our  analysis  we  are  going  to  take  into  account  the  third  order 
terms  of  equation  (52): 

/i 

T-'^G^^\Tz)  = 

h 

.  /4  . 

Where, 

fi  =  niiGf\Tz)+nnG2i3)iTz)  +  mzG^i\Tz)  (91) 

/2  =  n2iGf\Tz)  +  n22G2{S){Tz)  +  n23G^i\Tz)  (92) 

/3  =  n3,Gf\Tz)+ns2G2{S){Tz)  +  n33G^^\Tz)  (93) 
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h  =  n4iG^?{Tz)+n42G2i3)iTz)  +  n4iGf\T:z)  (94) 

If  we  substitute  equations  (68),  (69),  (70),  (71)  and  the  expressions  for  Gi, 
G2,  G3  into  equations  (91)  throuh  (94),  we  get: 

/l  —  1\2z\z2  +  l\zZ\Z2  +  I14Z2) 

+  ni2{l2iz\  +  ^22-2^1  •2^2  +  h3ZlZ2  +  I24Z2) 

+  ^13(^3l2^1  +  h2z‘\Z2  +  lzZZlZ2  +  Haz^)  (95) 

/2  =  ^21(^11-21  +  ^12-21 -22  +  I1ZZ1Z2  +  I14Z2) 

+  ^22(^21-21  +  I22Z1Z2  +  I2ZZ1Z2  +  /24-22) 

+  f^zzihiZi  +  lz2ZiZ2  +  I3ZZ1Z2  +  IZ4Z2)  (96) 

fz  —  ‘f^ZlihlZi  +  h2ZiZ2  +  I1ZZ1Z2  +  I14Z2) 

+  nz2{hlZi  +  l22z\z2  +  I2ZZ1Z2  +  I24Z2) 

+  ’233(^31-21  +  h2ZiZ2  +  hzziz^  +  Z342I)  (97) 

Ia  =  '>^4l{hlZi  +  Ii2Z^Z2  +  I1ZZ1Z2  +  I14Z2) 

+  n42(^21-2i  +  h2ZiZ2  +  hzZlZ^  +  ^242!) 

+  '’^4z{hiz\  +  lz2ZiZ2  +  lzzZlZ2  +  ^3422)  (98) 

From  equation  (52)  and  from  the  system  matrix,  the  derivatives  of  the  first 
two  modified  coordinates,  ii  and  Z2  can  be  written  as  follows: 

-21  =  {ae)zi  —  (uq  +  uj'€)z2  +  FFi(zi,  Z2)  (99) 

-22  =  (u)o  F  c^'e)zi  +  a'ez2  F  FF2(zi,  Z2)  (199) 

where. 


FFi(2i,2:2)  =  di  +  /i  +  ei 


(191) 
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FF2{zi,  22)  =  <^2  +  /2  +  62 


(102) 


If  we  combine  equations  (101)  and  (102)  with  equations  (72)  through  (75), 
and  (95)  through  (98),  we  get: 

-^•^1  (•2^1) -2^2)  =  F  V12Z1Z2 ri2ZiZ2 ri4Z2 

+  Pu4  +  P12Z1Z2  +  pnzl  +  ei  (103) 

FF2{zi,  Z2)  =  r2i2i  +  r22zfz2  +  r23Zizl  +  r24zl 

+  P21Z1  +  P22Z1Z2  +  P23^2  +  (104) 

where  coefficients  rij  and  pij  are: 

^11  =  ’^11^11  +  13^31 

ri2  =  ni\li2  +  ni2^22  +  ^13^32 
ri3  =  niiZi3  +  ni2l2Z  +  ^13^33 
ri4  =  +  7112^24  +  ^13^34 

^21  =  7121/11  +  7122/21  +  7123/31 
7*22  =  7X21/12  +  7X22/22  +  7X23/32 
7’23  =  7121/13  +  7X22/23  +  7X23/33 
7'24  =  7X2i/i4  +  7X22/24  +  7X23/34 

or  generally, 

7'ij  TXj^/lj  ~i”  '^i2^2j  71^3/3^  X  Ij  2  j  1,  3,  4  (105) 

also, 

Pll  =  Tlli/is  +  7X12/25  +  7I13/35 
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Pl2  =  ^11^16  +  ^12^26  +  ^13^36 
Pl3  =  Wii/i7  +  ni2^27  +  ^13^37 
P21  =  7121^15  +  7122^25  +  7723^35 
i722  =  7721^16  +  7722^26  +  7723/36 
P23  =  772i/i7  +  7722/27  +  7723/37 

or  generally, 


Py  —  77ii/i/:  +  77i2/2fc  +  77i3/3)fc  7  =  1,2  ji  =  1,  2,  3  A:  =  J  +  4 


(106) 


The  coefficients  kj  i  =  1, 2,  Z  ;  =  1, 2, 3, 4  are  as  follows: 


/ii  = 


<7ll 


^[“^(£'o777n  +  3Eimlim2i  +  3E2miim2i  +  Esmh) 


1 


Ol2rC'n. 


T>  ^  67 
1 


-(Eimfi  +  3£'2777iim2i  +  3jE'3miim|i  +  £'477721) 


+  -  ZBB)m 

<7ll  .Cd, 


3 

41 


(107) 


/12  —  (3£o777ii7r7i2  +  3£i(r77iirr722  +  2r77iir77i27772l) 


+  3£2(777i27772i  +  27772177722777 n)  +  3Esml^m22)  +  ^{W  -  B)3Tnl-^m42] 

^12  C  D 

~  ^-®2((777iir7722  +  2r77ii777 1277721) 

4-  3£3((777i27772i  +  27n2l77722777ii)  +  3£477l2i77722) 

+  -  Xs£)37774i77742]  +  “  2b5)37774i77742  (108) 


/l3  = 


Qj  ^  r  ,  2  /  2 

L  (3£^0^11^12  3£/i(??2 22^^21  +  2??2iimi2^22) 


D  ^  67 

+  3£2(77722777ii  +  27772177722777 12)  +  3£3r7721 77X22) 

+  ^{W  -  £) 37774177742] 
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~  ~^[~^{^Eimxim\2  +  ZE2{m\2m2\  +  2miimi2m22) 

+  ZEz{m%2fTiii  +  2m2im22Wi2)  +  ZE^m2im\^ 

+  ^{^gW  -  XBB)Zm4iml2]  +  -  ZBB)Sm4iml2 

^14  —  (-^0^12  +  ZEimi-^m2i  +  3£'2^  12^22  ^^3^22) 

+  \(W-B)ml,\ 

^12  ^  D 

-  +  3E2mfjm2i  +  3^3772120122  +  £^477132) 

+  g(37GW^  -  7Cs£)772t2]  +  ^(2g1^  -  2s£)77242 

^21  = - J"["^“(£'0777n  +  3£l7?lii77t2i  +  3£277ill777|i  +  £^377121) 

+  ^(W-B)mi] 


(109) 


(110) 


022rC'Dj 

-d^IT 


+  3E2mfim2i  +  SE’smum^i  +  £^477121) 


+  g((CGVF  -  a;sB)m4i]  +  -  ZBB)m\^ 

0-21  C!  £) 

I22  —  ^[■■j^^(3.£'0^ii^i2  +  3£'i(m^2^22  +  2mii7ni2W2i) 


(111) 


+  ZE2{mi2'm\i  +  2m2im22mii)  +  ZEzm\^m22)  +  -{W  —  5)377744777.42] 

~  3“  ZE2{^{rn\ivri22  +  2mii777i2^2i) 

+  ‘^Bz{{m\2m‘2i  +  2777,2i77722771ii)  +  ZE4m\-^m22) 

+  ^{xgW  -  XBB)3mlxm42]  +  ^(^gW  -  ZBB)Zml^m42  (112) 


Q21  ,Cd^ 
’  5  ^  67 


(35o7?7ii7?742  +  35i  (77742"i21  +  2777ii777 1277722) 


+  352(77712^11  +  2777217772277712)  +  3537772177722) 

+  ^{W  -  5)377741777^2] 

~  ”^[”0^(^-®1^11^12  352(7771277721  +  2777117771277722) 
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+  3£^3(^22’^ii  +  2m2im22'mi2)  +  SE4m2im22) 


+  -  XBB)Sm4iml2]  +  -  ZBB)Zm4im\2 

^21  ^  D 

+  ‘iEim\-jm2i  +  3£2^12^22  + 


I24  = 


D  ‘  67 


+  -{W-B)Tn 


42J 


-  +  ^E2ml^m2i  +  3Ezmi2ml2  + 


+  -xi^G^  -  3:BB)ml2]  +  ^{zgW  -  ZBB)m% 


6D 


''42 


^31  =  — 3£'imfim2i  +  3E2miiml-^  +  Ezmh) 

+  l{W-B)ml] 

^32  C^D 

-  ^  3E2ml-^m2i  +  SEsmuml-j^  +  £'477121) 


1 


+  -(xgPF  -  XB£)m4i]  +  -^{zcW  -  ZBB)m 
_  031 


61) 
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^32  =  - 


D  ‘  67 


(3£omiimi2  +  3£i(miim22  +  2mu'mi2m2i) 


(113) 


(114) 


(115) 


+  3£2(mi2m2i  +  2m2im22'mii)  +  ZEzm\im22)  +  q(^  ~  B)3m^im42] 

-  -^[-^(3£imiimi2  +  3£2((miim22  +  2miimi2m2i) 

+  3Ez{{mi2mli  +  2m2im22mii)  +  3E4m\im22) 


1 


+  q{xgW  -  XBB)3ml^m42\  +  ^i^G^  -  2b£) 37714177142 
hz  =  —^[-^(3£omii771i2  +  3£i(771i277l21  +  277111771 1271122) 
+  3£2(77l22771ii  4-  2711217712277112)  +  3Ezm2im%^ 


(116) 


+  g(^4^  -  £)3r7i4i77iy 


Q32  fCp^ 
D  ^  67 


(3£ir71ii771i2  +  3£2(771i277l21  +  2771iir71i277l22) 


+  3£3(77l22771ll  +  277l2l77l22771i2)  +  3£47?l2l77l22) 
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(117) 


+  ^{xaW  -  iBS)3m4im|2)  +  -  XBB)3m4,mh 

^31  ^jD 

I34  =  ~"^I“g~(^'0*^i2  +  3£^imf2m2i  +  3£2^i2^22  +  ^3^32) 

+  i(M' -B)mi] 

“  "^[“^^(^1^12  +  3£^2?^11^21  +  3Esmi27n22  +  £^4/7122) 

+  q(^g^  -  2:55)77142]  +  -  ZBB)m\2  (118) 

The  next  step  is  to  introduce  polar  coordinates  in  the  form: 

21  =  5  cos  ^  (119) 

22  =  5  sin  0  (120) 

We  use  polar  coordinates,  because  it  is  easier  in  this  way  to  investigate  the 
existence  of  limit  cycles. 

Substituting  equations  (119),  (120),  into  equations  (99),  (100),  we  get: 

R  cos  6  —  R{sixL  0)6  =  (a;o  +  ai^e)5cos0  +  o^ei?sin0  + 

Pi{e)R^  +  Q^{e)R^  +  ex  (121) 

R  sm  0  +  R{cos  9)0  =  (wq  +  a;'e)5  cos  6^  +  a'e5  sin  0 + 

P2{e)R^  +  Q2{e)R^  +  62  (122) 

where, 

Pi(^)  =  rii  cos^  6*  +  ri2  cos^  0  sin  0  +  ri3  cos  0  sin^  ^  +  ri4  sin^  6'  (123) 

Qi{^)  =  Pii  cos^0 +  Pi2cos^sin0 +Pi3sin^0  (124) 

P2{d)  =  r2i  cos^  6*  +  r22  cos^  6*  sin  0  +  r23  cos  6>  sin^  0  +  r24  sin^  6*  (125) 

Q2{0)  =  P21  cos^  0  +  P22  cos  6>  sin  ^  +  P23  sin^  ^  (126) 
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If  we  multiply  equation  (121)  by  cos^  and  equation  (122)  by  sin6>,  and  add 
the  two  resulting  equations,  we  get; 


R  —  a  eR  +  P{d)R^  +  Q(d)R?‘  +  (ei  cos  6  +  e2  sin  6)  (127) 


where, 


P{9)  =  Pi{9)cos6  +  P2{9)sm9  (128) 

<3(^)  =  <5i(^)cos^  +  (52(^)sin6>  (129) 

Equation  (127)  contains  one  variable  that  varies  slowly  in  time  {R)  and  a  fast 
variable  {9). 

If  we  average  this  equation  over  one  complete  cycle  in  9,  from  0  to  27r, 
equation  (127)  becomes: 


R  =  a'eR  +  KR^  +  LRj^  +  M  (130) 


where. 


1  /•2’r 

^ 

—  g(3rii -f  ri3  +  r22  +  3r24  (131) 

^  (132) 

1 

M  =  —  /  d  cos  9  d9 
2tt  Jo 

r2TT 

+  —  62  sin  9  d9 

27r  Jo 

=  |^[cos0]2’i 

=  0  (133) 
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Finally  equation  (135)  becomes: 


R  =  aeR  +  KR^  (134) 

F.  LIMIT  CYCLE  ANALYSIS 

At  steady  state  R  =  0,and  equation  (134)  becomes: 

Q^R{a'e  +  KR?)  (135) 

Equation  (135)  has  two  solutions.  The  first  solution  is  i?  =  0.  This  is  the 
trivial  solution  and  it  does  not  give  us  much  information. The  second  solution 
is: 

(136) 

This  solution  gives  us  a  limit  cycle  of  constant  amplitude  R  in  the  zi, 
cartesian  coordinate  system.This  limit  cycle  exists  if  the  quantity  inside  the 
square  root  is  possitive,  or 

—a'e 

—  >  0  (137) 

Condition  (137)  is  necessary  for  the  amplitude  of  the  limit  cycle,  i?,  to  be  a 
real  number. 

In  our  case  a'  is  always  negative,  because  for  constant  rc^,  as  we  decrease 
e  (Figure  1),  the  real  part  of  the  critical  pair  of  eigenvalues  increases,  due  to 
further  loss  of  stability.  In  other  words  we  can  say  that: 

a'  <  0  (138) 
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From  conditions  (137),  (138)  we  see  that  the  existence  of  the  limit  cycles 
depends  on  the  value  of  parameter  K.  We  can  see  that: 

•  UK  <0,  periodic  solutions  exist  for  e  <  0  or  zq  —  zg^  <  Q  oi  zg  <  zg^, 
and  they  are  stable. 

•  UK  >0,  periodic  solutions  exist  for  e  >  0  or  zg  -  zg^  >  0  or  zg  >  zg^, 
and  they  are  unstable. 

The  characteristic  root  of  equation  (134)  in  the  vicinity  of  (136)  is: 

0  =  -2a'e  (139) 

The  sign  of  this  characteristic  root,  assigns  the  stability  of  the  periodic  solu¬ 
tions. 


G.  RESULTS  AND  DISCUSSION 

Typical  results  in  terms  of  the  nonlinear  stability  coefficient  K  are  pre¬ 
sented  in  Figures  3  through  6.  The  stability  coefficient  K  is  shown  in  its 
normalized  form  as  A  •  7  (Papadimitriou,  1994).  Figure  3  shows  K  versus 
the  LCG/LCB  separation  distance  xg  (in  ft)  for  a  given  vehicle  speed  and  for 
different  values  of  the  drag  coefficient.  It  can  be  seen  that  K  is  everywhere 
negative,  which  means  that  all  bifurcations  to  periodic  solutions  are  super¬ 
critical.  Higher  values  of  the  drag  coefficient  result  in  stronger  supercritical 
bifurcations  which  means  that  the  corresponding  limit  cycle  amplitudes  will 
be  smaller. 
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Figure  4  shows  K  versus  xq  for  a  given  drag  coefficient  and  for  different 
forward  speeds.  It  can  be  seen  that  the  bifurcations  are  supercritical  with  the 
possible  exception  of  large  speeds  and  small  xq.  This  means  that  it  is  possible 
for  a  properly  trimmed  vehicle  at  relatively  high  speeds  to  experience  an  oscil¬ 
latory  behavior  even  before  stability  is  lost.  This  demonstrates  a  destabilizing 
effect  which  could  not  have  been  predicted  by  linear  techniques. 

Figures  5  and  6  shows  numerical  simulation  results  which  confirm  the  theo¬ 
retical  predictions.  Both  figures  correspond  to  a  supercritical  bifurcation  case 
and  they  show  a  continuous  increase  in  limit  cycle  amplitudes  as  the  bifurca¬ 
tion  point  is  crossed. 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 


This  thesis  presented  a  comprehensive  nonlinear  study  of  straight  line  sta¬ 
bility  of  motion  of  submersibles  in  coupled  sway /yaw /roll  motions  under  open 
loop  conditions.  Primary  loss  of  stability  was  shown  to  occur  in  the  form  of 
Hopf  bifurcations  to  periodic  solutions.  This  loss  of  stability  is  characteristic 
of  the  coupling  of  roll  into  sway  and  yaw  and  cannot  be  predicted  by  consid¬ 
ering  the  uncoupled  motions.  The  critical  point  where  instability  occurs  was 
computed  in  terms  of  vehicle  metacentric  height,  longitudinal  separation  of 
the  centers  of  buoyancy  and  gravity,  and  the  forward  speed.  Analysis  of  the 
periodic  solutions  that  resulted  from  the  Hopf  bifurcations  was  accomplished 
through  Taylor  expansions,  up  to  third  order,  of  the  equations  of  motion.  A 
consistent  approximation,  utilizing  the  generalized  gradient,  was  used  to  study 
the  non-analytic  quadratic  cross  flow  integral  drag  terms.  The  results  indi¬ 
cated  that  loss  of  stability  occurs  always  in  the  form  of  supercritical  Hopf 
bifurcations  with  stable  limit  cycles.  It  was  shown  that  this  is  mainly  due  to 
the  stabilizing  effect  of  the  drag  forces  at  high  angles  of  attack.  Subcritical 
bifurcations,  however,  with  considerably  higher  limit  cycle  amplitudes  may 
develop  for  sufficiently  high  forward  speeds  and  small  LCG/LCB  separations. 
Simulation  studies  of  these  subcritical  bifurcations  along  with  the  effects  of 
vertical  plane  coupling  constitute  recommendations  for  further  research. 
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APPENDIX 


The  following  is  a  list  and  description  of  the  computer  programs  used  in  this 
thesis.  The  programs  are  written  in  FORTRAN  or  MATLAB.  Complete  print¬ 
outs  of  the  programs  follow  after  the  list. 

•  STABILITY.M 

MATLAB  program  for  performing  linear  stability  analysis. 

•  SIM.M 

MATLAB  program  and  functions  for  numerical  simulation. 

•  HOPF.FOR 

FORTRAN  program  for  calculating  the  nonlinear  stability  analysis  coef¬ 
ficient  K.  It  requires  data  from  STABILITY.M  and  standard  numerical 
linear  algebra  subroutines. 
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•/.  STABILITY. M 

7. 

7.  LOSS  OF  STABILITY 
a=l 

W=12000; 

IXX=1760;  IYY=9450; 

IZZ=10700:  IXZ=0;  IXY=0; 

IYZ=0;  L=17.425;  RH0=1.94; 

G=32.2;  U=6.5;  M=W/G;  B=W; 

ND1=0.5*RH0*L~2; 

7.  DEFINE  HYDRODYNAMIC  COEFFICIENTS 

YPDOT=l .  270e-04*NDl*L-'2 ; 

YVD0T=-5 . 550e-02*NDl*L ; 

YRDOT=l .  240e-03*NDl*L''2 ; 
YP=3.055e-03*NDl*L; 
YV=-9.310e-02*NDl; 
YR=-5.940e-02*NDl*L; 

NPD0T=-3 . 370e-05*NDl*L~3 ; 
NVD0T=1.240e-03*NDl*L*2; 

NRD0T=-3 . 400e-03*NDl*L~3 ; 

NP=-8 . 405e-04>*=NDl*L~2 : 

NV=-1 .484e-02*NDl*L; 

NR=-1 .640e-02*NDl*L''2; 

KPD0T=-1.01e-03*NDl*L'3; 

KVDOT=l . 27e-04*NDl*L~2 ; 

KRD0T=-3 . 37e-05*NDl*L‘3 ; 

KP=-1 . 10e-02*NDl*L‘2; 
KV=3.055e-034=NDl*L; 

KR=-8 . 41e-04’*'NDl*L~2 ; 

flag=0; 

for  XG=0:0.01:2, 
f lag=f lag+1 ; 
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xg(flag)=XG; 

a=IXX-KPDOT;  b=KP*U;  e=KV=t:U; 
f=KRDOT;  i=YP*U;  j=M-YVDOT;  k=YV*U; 

1=XG*M-YRD0T;  ni=U*(YR-M);  o=NPDOT; 
p=NP*U;  q=-XG*W;  r=XG*M-NVDOT; 
w=U=t=  (NR-XG*M)  ;  x=NV*U;  u=IZZ-NRDOT: 

al=-u*M~2; 

a2=-u*M*YPD0T-u*M*KVD0T+M*o*l+f*r=t'M; 

a3=a* j  *u-a*l*r-u*KVDOT*YPDOT+KVDOT*o*l+f *r*YPDOT-f *0* j ; 

bl=(w*(M“2))+(r*U*(M~2)) ; 

b2=-M*e*u-M*u*i+w*M=i'YPD0T+w*KVD0T*M-o*m*M+p*l*M-f*x*M+. 

r*M*U*YPDOT-o*j*M*U+r*KR*U*M; 

b3=-e*u*YPD0T+e*o*l-u*i*KVD0T+w*KVD0T*YPD0T-KVD0T=t=o=t:m+KVD0T=i=p*l 
a*  j  *w-a*k*u+a*l*x+a*r=i'in-b*  j  *u+b*l^'r+f  =('r*i-f  *x*YPD0T+ . 
o*k*f-f *p*j+r*KR*U*YPDOT; 

cl=-x*(M~2)=t'U; 

c2=r*i*M*U-x*M*U*YPD0T-x*KR*U*M+o*k*M*U-p* j *M*U+ j *u*W-l*r*W- 
q*l*M+w*i*M-m*p*M+e*w*M; 

c3=a*k*w-a*m*x+b* j  *w+b*k*u-b*l*x-b*r*m+r*i*KR*U-x*KR*U*YPDOT+ 
o=t=k*KR*U-p*j*KR*U-q*l*KVDOT+w*i*KVDOT-in*p*KVDOT-e*u*i+e*w*YPDOT 

e+o^m+e+p’t'l+f  =t:q*  j  -f  =t:x=t'i+f  ^p^k ; 

d2=q*j*M*U-x*i*M*U+p*k*M=i'U+q*ni*M-j*w*W-k*u*W+l*x*W+r*m*W; 
d3=q*j*KR*U-x+i*KR*U+p*k*KR*U-f*q*k+q*m=i'KVD0T-e*q*l+e*w*i- . 
e*m*p-b*k*w+b*m*x ; 

e2=k*w*W-m*x*W-q*k=t'M=i'U ; 
e3=e*q*m-q*k*KR*U ; 

f5=(-e2=t=(bl-2))+bl*cl*d2; 

f4=((b2*cl+bl*c2)*d2)+bl*cl*d3-al*(d2''2)-e3*(bl~2)-2*e2*bl*b2; 
f3=-a2=t=(d2~2)-2*d3*d2*al-2*e3*bl=t=b2-e2*^((b2~2)+2*bl*b3)  +  . .  . 
d2* (b3*cl+b2*c2+bl*c3)+d3* (b2*cl+bl*c2) ; 
f2=-e3*((b2''2)+2*bl*b3)-2*e2*b2*b3+d2*(b3*c2+b2*c3)+. . . 
d3=t=(b3=t=cl+b2*c2+bl*c3)-a3*(d2*2)-2>t:d3*d2*a2-al*(d3*2) ; 
f I=b3*c3*d2+d3*(b3*c2+b2*c3)-2*e3*b2*b3-e2*(b3“2)-. . . 
2*d3*d2*a3-a2*(d3'2) ; 
f 0=b3*c3*d3-a3* (d3~2) -e3* (b3~2) ; 


55 


coef=[f5  f4  f3  f2  fl  fO] ; 

ZG=roots(coef ) ; 

tot(flag)=ZG(5,l) ; 
end 

plot (xg.tot) .grid; 

title (’ZG  versus  XG  plot  in  the  point  of  loss  of  stability’) 
xlabeK’XG  in  ft’); 
ylabeK’ZG  in  ft’) ; 
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7.  NON  LINEAR  SIMULATION  PROGRAM 
t0=0; 

tfinal=500; 
q=(l/180)*pi; 
y0=[0  0  0  q]  ; 

[t , y] =ode45 ( ' vdpo2 ’ , tO , tf inal , yO) ; 

figure (1) ,plot(t(: , 1) ,57. 29578*y( ; ,4) ) ,grid 
titleC’plot  of  fi  with  time’); 
ylabel(’fi’) ,xlabel(’time’) 


57 


%  NON  LINEAR  SIMULATION  PROGRAM 

f miction  5rprime=vdpo2(t,y) 

W=12000; 

IXX=1760;  IYY=9450; 

IZZ=10700;  IXZ=0;  IXY=0; 

IYZ=0;  L=17.425;  RH0=1.94;  ZB=0; 

G=32.2;  U=5;  M=W/G;  B=W;  XB=0; 

ZG=0.05;  CD=0.5;  YG=0; 

XG=1;  YB=0; 

ND1=0.5*RH0*L'2; 

%  DEFINE  HYDRODYNAMIC  COEFFICIENTS 

YPDOT=l .  270e-04*NDl>^L‘'2 ; 

YVD0T=-5 . 550e-02*NDl*L ; 

YRD0T=1.240e-03*NDl*L‘2; 

YP=3.055e-03*NDl*L; 

YV=-9.310e-02*NDl; 

YR=-5.940e-02*NDl*L; 

NPD0T=-3.370e-05*NDl*L~3; 

NVD0T=1 . 240e-03*NDl*L''2 ; 

NRDOT=-3 . 400e-03*NDl*L~3 ; 

NP=-8 . 405e-04*NDl*L‘2 ; 

NV=-1.484e-02*NDl*L; 

NR=-1.640e-02*NDl*L~2; 

KPD0T=-1.01e-03*NDl*L'3; 

KVDOT=l .  27e-04*NDl=»=L'2 ; 

KRD0T=-3 . 37e-05*NDl*L*3 ; 

KP=-1 . 10e-02*NDl*L‘2; 

KV=3.055e-03*NDl*L; 

KR=-8.41e-04*NDl*L''2; 

D=((M-YVDOT)*(IZZ-NRDOT)*(IXX-KPDOT)) . . . 
-((M-YVDOT)*(IXZ-KRDOT)*(-IXZ-NPDOT)) . . . 

-  (  (M*XG-NVDOT)  =*=  (M*XG-YRDOT)  *  (IXX-KPDOT) )  . . 
+ ( (M*XG-NVDOT) * (IXZ-KRDOT) * (-M*ZG-YPDOT) ) . 
+ ( (-M*ZG-KVDOT) * (M+XG-YRDOT) * (-IXZ-NPDOT) ) 
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- ( (-M*ZG-KVDOT) * (IZZ-NRDOT) * (-M*ZG-YPDOT) ) ; 

All= ( (IZZ-NRDOT) * (IXX-KPDOT) ) - ( (IXZ-KRDOT) * (-IXZ-NPDOT) ) ; 

A12= ( (-M*XG+YRDOT) * ( IXX-KPDOT) ) + ( (IXZ-KRDOT) * (-M+ZG-YPDOT) ) ; 
A13=  (  (M*XG-YRDOT)  *  (-IXZ-NPDOT)  )  -  (  (IZZ-NRDOT)  *  (-M*ZG-YPDOT) ) 
A21=  (  (-M*XG+NVDOT)  *  (IXX-KPDOT) )  +  (  (-M=t=ZG-KVDOT)  *  (-IXZ-NPDOT) )  ; 
A22= ( (M-YVDOT) * (IXX-KPDOT) ) - ( (-M*ZG-KVDOT) * (-M+ZG-YPDOT) ) ; 

A23= ( (-M+YVDOT) * (-IXZ-NPDOT) ) + ( (M*XG-NVDOT) * (-M*ZG-YPDOT) ) ; 
A31= ( (M*XG-NVDOT) * (IXZ-KRDOT) ) - ( (-M*ZG-KVDOT) * (IZZ-NRDOT) ) ; 
A32= ( ( -M+YVDOT ) * ( IXZ-KRDOT ) ) + ( ( -M*ZG-KVDOT) * (M*XG-YRDOT) ) ; 

A33= ( (M-YVDOT) * (IZZ-NRDOT) ) - ( (M*XG-NVDOT) * (M*XG-YRDOT) ) ; 

7.  EVALUATE  TRANSFORMATION  MATRIX  OF  EIGENVECTORS 

% 

F (1 , 1) = (A11*YV*U+A12*NV*U+A13*KV*U) /D ; 

F(1,2)=(A11* (YR*U-M*U) +A12* (-M*XG*U+NR*U) +A13* (M*ZG*U+KR*U) ) /D ; 
F (1 , 3) = (A11*YP*U+A12*NP*U+A13*KP*U) /D ; 

F(l,4)=0; 

F (2 , 1 ) = (A21 *YV*U+A22*NV*U+A23*KV*U) /D ; 

F  (2 , 2)  =  (A21*  (YR*U-M*U)  +A22*  (-M*XG*U+NR*U)  +A23*  (M>t'ZG*U+KR*U) )  /D ; 
F  (2 , 3)  =  (A21*YP*U+A22>i'NP*U+A23*KP*U)  /D ; 

F(2,4)=0; 

F (3 , 1) = (A31*YV*U+A32*NV*U+A33*KV*U) /D ; 

F  (3 , 2)  =  (  A31*  (  YR*U-M=*=U)  +A32*  (-M*XG*U+NR=t=U)  +A33*  (M*ZG*U+KR+U)  )  /D ; 
F (3 , 3) = (A31*YP*U+A32*NP*U+A33*KP*U) /D ; 

F(3,4)=0; 

F(4,l)=0; 

F(4,2)=0; 

F(4,3)=l: 

F(4,4)=0: 

7.  CALCULATION  OF  THE  INTEGRATION  TERMS  Ei 


7.  DEFINE  THE  LENGTH  X  AND  HEIGHT  H  TERMS  FOR  THE  INTEGRATION 
7.  ALL  IN  FEET. 

XL(1)=-105.9/12; 

XL(2)=-104.3/12; 

XL(3)=-99.3/12; 

XL(4)=-94.3/12; 

XL(5)=-87.3/12; 

XL(6)=-76.8/12; 
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XL(7)=-66.3/12; 

XL(8)=-55.8/12; 

XL(9)=72.7/12: 

XL(10)=79.2/12; 

XL(11)=83.2/12; 

XL(12)=87.2/12; 

XL(13)=91.2/12; 

XL(14)=95.2/12; 

XL(15)=99.2/12; 

XL(16)=101.2/12; 

XL(17)=102.1/12; 

XL(18)=103.2/12; 


HT(1)=0: 
HT(2)=2. 28/12; 
HT(3)=8. 24/12; 
HT(4)=13. 96/12; 
HT(5)=19. 76/12; 
HT(6)=25.1/12; 
HT(7)=29. 36/12; 
HT(8)=31. 85/12; 
HT(9)=31. 85/12; 
HT(10)=30. 00/12 
HT(11)=27. 84/12 
HT(12)=25. 12/12 
HT(13)=21. 44/12 
HT(14)=17. 12/12 
HT(15)=12.0/12; 
HT(16)=9. 12/12; 
HT(17)=6. 72/12; 
HT(18)=0; 


for  i=l:18, 

VECl(i)=HT(i)*(y(l)+XL(i)*y(2))=f=(abs(y(l)+XL(i)*y(2))); 
VEC2(i)=HT(i)>^(y(l)+XL(i)*y(2))=t=(abs(y(l)+XL(i)*y(2)))=t'XL(i)  ; 
end 


0UT=0; 
for  j=l:17, 

0UT1=0 . 5* (VECl ( j ) +VEC1 ( j +1) ) * (XL ( j +1 ) -XL ( j ) ) ; 

0UT=0UT+0UT1 ; 

end 
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IV=CD*OUT: 


0UT=0; 
for  j=l:17, 

0UT2=0 . 5* (VEC2 ( j ) +VEC2 ( j +1 ) ) * (XL ( j  +1 ) -XL ( j ) ) ; 

0UT=0UT+0UT2; 

end 

IR=CD*0UT; 


yprime=[F(l,l)*y(l)+F(l,2)*y(2)+F(l,3)*y(3)+F(l,4)*y(4)+. . . 
(All/D)*(YG*(y(3)*2+y(2)~2)-IV+(W-B)*sin(y(4)))+. . . 
(A12/D)*(IXY*y(3)~2+IYZ+y(3)*y(2)+YG+y(l)*y(2)-IR+(XG*W-XB*B)*sin(y(4)))+ 
(A13/D)*(-IXY*y(3)*y(2)-IYZ*y(2)'‘2-M*YG*y(l)*y(3)+(YG*W-YB*B)*cos(y(4)) 
-(ZG*W-ZB*B)*sin(y(4))); . . . 

F(2,l)*y(l)+F(2,2)*y(2)+F(2,3)*y(3)+F(2,4)*y(4)+. . . 
(A21/D)*(YG*(y(3)-'2+y(2)'2)-IV+(W-B)*sin(y(4)))  +  .  . . 
(A22/D)*(IXY=t=y(3)'‘2+IYZ=('y(3)*y(2)+YG*y(l)*y(2)-IR+(XG*W-XB*B)*sin(y(4)))  + 
(A23/D)*(-IXY*y(3)*y(2)-IYZ=i'y(2)~2-M*YG*y(l)*y(3)  +  (YG*W-YB*B)*cos(y(4)) .  . 
-(ZG*W-ZB*B)*sin(y(4)))  ; .  .  . 

F(3,l)*y(l)+F(3,2)*y(2)+F(3,3)*y(3)+F(3,4)*y(4)  +  . . . 
(A31/D)*(YG*(y(3)''2+y(2)~2)-IV+(W-B)*sin(y(4)))  +  .  .  . 
(A32/D)*(IXY*y(3)‘2+IYZ*y(3)*y(2)+YG*y(l)*y(2)-IR+(XG*W-XB*B)*sin(y(4)))+ 
(A33/D)*(-IXY*y(3)=t=y(2)-IYZ*y(2)“2-M*YG*y(l)*y(3)  +  (YG*W-YB=t=B)*cos(y(4)) . . 
-(ZG*W-ZB*B)*sin(y(4)))  ; .  .  . 
l*y(3)] ; 
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HOPF.FOR 


C 
C 

C  HOPF  BIFURCATIONS  PROGRAM 
C 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

REAL*8  L , lYY , M , YPDOT , YVDOT , NDl , YRDOT 
REAL*8  YP , YV , YR, NPDOT , NVDOT , NRDOT , NP , NV , NR 
REAL*8  GAMA,U,KVDOT,KRDOT,KPDOT,D 
REAL*8  E0,E1,E2,E3,E4,XG,ZG,KR,KP,KV,XG1,ZG1 
C 

REALMS  M11,M12,M13,M14,M21,M22,M23 
REAL*8  M24 , M31 , M32 , M33 , M34 , M41 , M42 , M43 , M44 
REAL*8  N11,N12,N13,N14,N21,N22,M23,N24 
REAL*8  N3 1 , N32 , N33 , N34 , N4 1 , N42 , N43 , N44 
REAL*8  L11,L12,L13,L14,L21,L22,L23.L24,L31 
REAL*8  L32 , L33 , L34 , L15 , L16 , LIT , L25 , L26 , L27 , L35 
REAL*8  L36 , L37 , LI A , L2A , L3A , L4A , L5A , L6A , L7A , L8A , L9A 
REAL*8  L10A,L11A,L12A,L1,L2,L3,L4,L5,L6,L7 
REAL*8  L8,L9,R11,R12,R13,R14,R21,R22,R23,R24 
REAL*8  P11,P12,P13,P21,P22,P23 
C 

DIMENSION  F (4 , 4) , T (4 , 4) , TINV (4 , 4) , FVl (4) , IVl (4) , YYY (4,4) 
DIMENSION  WR(4) ,WI(4) ,TSAVE(4,4) ,TLUD(4,4) , IVLUD(4) ,SVLUD(4) 
DIMENSION  ASAVE(4,4) ,A2(4,4) ,XL(18) ,HT(18) ,ZGG(197) ,FF(4,4) 
DIMENSION  VEC0(18) ,VEC1(18) ,VEC2(18) ,VEC3(18) ,VEC4(18) ,XGG(197) 
C 

INTEGER  I,J,K 
C 

OPEN  ( 20 , FILE= ’ HOPF .RES’, STATUS= ’ OLD ’ ) 

OPEN  (2 1 , FILE= ’ DATA . DAT ’ , STATUS= ’ OLD ’ ) 

OPEN  ( 23 , FILE= ’ HOPF 1 . RES ’ , STATUS= ’ OLD ’ ) 

OPEN  (25,FILE=’H0PF2.RES’ ,STATUS=’OLD’) 

C 


WEIGHT=12000 . 0 


IXX 

=1760.0 

lYY 

=9450.0 

IZZ 

=10700.0 

IXZ 

=0.0 

IXY 

=0.0 

lYZ 

=0.0 

L 

=17.425 

RHO 

=1.94 
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WRITE  (*,*)  ’  ENTER  CD’ 

READ  (*,*)  CD 
G  =32.2 

XB  =0.0 

ZB  =0.0 

YG  =0.0 

YB  =0.0 

YDELTAR=0.0 
DELTAR=0.0 
NDELTAR=0.0 
NPR0P=0 . 0 
M=  WEIGHT/G 
B=  WEIGHT 
W=WEIGHT 
U=8.0 

ND1=0.5*RH0*L**2 

C 

C  DEFINE  HYDRODYNAMIC  COEFFICIENTS 
C 

YPD0T=1 . 270E-04*ND1*L**2 
YVD0T=-5 . 550E-02*ND1*L 
YRD0T=1 . 240E-03*ND1*L**2 
YP=3.055E-03*ND1*L 
YV=-9.310E-02*ND1 
YR=-5.940E-02*ND1*L 
C 

NPD0T=-3 . 370E-05*ND1*L**3 
NVD0T=1 . 240E-03*NDl*L*=i'2 
NRD0T=-3 .400E-03*ND1*L**3 
NP=-8 .405E-04=t=NDl*L**2 
NV=-1 .484E-02*ND1*L 
NR=-1 . 640E-02=t=NDl*L=*=*2 
C 

KPD0T=-1 . 01E-03*ND1*L**3 
KVD0T=1 . 27E-04+NDl*L=t=*2 
KRD0T=-3 . 37E-05*ND1*L**3 
KP=-1 . 10E-02*ND1*L**2 
KV=3.055E-03*ND1*L 
KR=-8 . 41E-04*ND1*L**2 
C 

C  DEFINE  THE  LENGTH  X  AND  HEIGHT  H  TERMS  FOR 
C  THE  INTEGRATION,  ALL  IN  FEET. 
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c 


XL(  1)=-105. 9/12.0 
XL(  2)=-104. 3/12.0 
XL(  3)=-99. 3/12.0 
XL(  4)=-94. 3/12.0 
XL(  5)=-87. 3/12.0 
XL(  6)=-76. 8/12.0 
XL(  7)=-66. 3/12.0 
XL(  8)=-55. 8/12.0 
XL(  9)=72. 7/12.0 
XL(10)=79. 2/12.0 
XL(11)=83. 2/12.0 
XL(12)=87. 2/12.0 
XL(13)=91. 2/12.0 
XL(14)=95. 2/12.0 
XL(15)=99. 2/12.0 
XL(16)=101. 2/12.0 
XL(17)=102. 1/12.0 
XL(18)=103. 2/12.0 
C 

HT(  1)=  0.000 
HT(  2)=  2.28/12.0 
HT(  3)=  8.24/12.0 
HT(  4)=  13.96/12.0 
HT(  5)=  19.76/12.0 
HT(  6)=  25.1/12.0 
HT(  7)=  29.36/12.0 
HT(  8)=  31.85/12.0 
HT(  9)=  31.85/12.0 
HT(10)=  30.00/12.0 
HT(11)=  27.84/12.0 
HT(12)=  25.12/12.0 
HT(13)=  21.44/12.0 
HT(14)=  17.12/12.0 
HT(15)=  12.0/12.0 
HT(16)=  9.12/12.0 
HT(17)=  6.72/12.0 
HT(18)=  0.00 
C 

DO  104  K  =  1,18 
VECO(K)=HT(K) 
VEC1(K)=XL(K)*HT(K) 
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VEC2 (K) =XL (K) *XL (K) *HT (K) 

VEC3 (K) =XL (K) *XL (K) *XL (K) *HT(K) 

VEC4 (K) =XL (K) *XL (K) *XL (K) *XL (K) *HT (K) 

104  CONTINUE 

CALL  TRAP(18,VEC0,XL,E0) 

CALL  TRAP(18,VEC1,XL,E1) 

CALL  TRAP(18,VEC2,XL,E2) 

CALL  TRAP(18,VEC3,XL,E3) 

CALL  TRAP(18,VEC4,XL,E4) 

C 

WRITE  (*,*)  '  ENTER  GAMA’ 

READ  (*,*)  GAMA 

C  READ  THE  CRITICAL  VALUES  FOR  XG  AND  ZG  FROM  FILE  DATA.DAT 

C 

XGG(1)=0.0 
ZGG(1)=0. 016358083 
DO  1  1=2,197 

READ  (21,*)XG,ZG 

XGG(I)=XG 

ZGG(I)=ZG 

C  DETERMINE  [F]  COEFFICIENTS 

C 

D=( (M-YVDOT) * (IZZ-NRDOT) * (IXX-KPDOT) ) 

&  - ( (M-YVDOT) * (IXZ-KRDOT) * (-IXZ-NPDOT) ) 

&  - ( (M*XG-NVDOT) * (M*XG-YRDOT) * (IXX-KPDOT) ) 

&  + ( (M*XG-NVDOT) * (IXZ-KRDOT) * (-M*ZG-YPDOT) ) 

&  + ( (-M*ZG-KVDOT) * (M*XG-YRDOT) * (-IXZ-NPDOT) ) 

&  - ( (-M*ZG-KVDOT) * (IZZ-NRDOT) * (-M+ZG-YPDOT) ) 

C 

A1 1= ( (IZZ-NRDOT) * (IXX-KPDOT) ) - ( (IXZ-KRDOT) * (-IXZ-NPDOT) ) 
A12=(  (-M=t=XG+YRDOT)*  (IXX-KPDOT)  )  +  (  (IXZ-KRDOT)  *  (-M*ZG-YPDOT)  ) 
A13= ( (M*XG-YRDOT) * (-IXZ-NPDOT) ) - ( (IZZ-NRDOT) ♦ (-M*ZG-YPDOT) ) 
A2 1= ( (-M^XG+NVDOT) * ( IXX-KPDOT) ) + ( (-M*ZG-KVDOT) * (-IXZ-NPDOT) ) 
A22=  (  (M-YVDOT)  =t=  ( IXX-KPDOT) )  -  (  (-M*ZG-KVDOT)  *  (-M*ZG-YPDOT)  ) 
A23=( (-M+YVDOT)*(-IXZ-NPDOT) )+((M*XG-NVDOT) *(-M*ZG-YPDOT) ) 
A31=( (M*XG-NVDOT) * (IXZ-KRDOT) ) - ( (-M*ZG-KVDOT) * (IZZ-NRDOT) ) 
A32=  (  (-M+YVDOT)  *  ( IXZ-KRDOT) )  +  (  (-M*ZG-KVDOT)  *  (M=(»XG-YRDOT)  ) 
A33=  (  (M-YVDOT)  *  (IZZ-NRDOT)  )  -  (  (M+XG-NVDOT)  *  (M=t=XG-YRDOT)  ) 

C 

F  ( 1 , 1)  =  (All*YV*U+A12=i'NV=*=U+A13*KV*U)  /D 
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F(1,2)  =  (A11*  (YR*U-M*U)  +A12*  (-M*XG*U+NR*U)  +A13*  (M*ZG=*=U+KR*U)  )  /D 
F(l,3)  =  (All*YP*U+A12*NP*U+A13*KP=t=U)  /D 
F(1,4)=(A11* (W-B) +A12* (XG*W-XB*B) +A13* (-ZG*W+ZB*B) ) /D 
F  (2 , 1)  =  (A21*YV*U+A22>t=NV*U+A23*KV*U)  /D 

F (2 , 2) = (A21* (YR*U-M*U) +A22* (-M*XG*U+NR*U) +A23* (M*ZG*U+KR*U) ) /D 
F  (2 , 3)  =  (A21*YP*U+A22*NP*U+A23=t=KP*U)  /D 
F (2 , 4) = (A21* (W-B) +A22* (XG*W-XB*B)+A23* (-ZG*W+ZB*B) ) /D 
F (3 , 1 ) = (A31*YV*U+A32*NV*U+A33*KV*U) /D 

F(3,2)=(A31*(YR*U-M*U)+A32*(-M*XG*U+NR*U)+A33*(M*ZG*U+KR*U))/D 
F (3 , 3) = (A31*YP*U+A32*NP*U+A33*KP*U) /D 
F  (3 ,4)  =  (A31*  (W-B)  +A32*  (XG*W-XB*B)+A33=i'  (-ZG*W+ZB*B) )  /D 
F(4,l)=0.0 
F(4,2)=0.0 
F(4,3)=1.0 
F(4,4)=0.0 
C 

C  EVALUATE  TRANSFORMATION  MATRIX  OF  EIGENVECTORS 

C 

DO  11  K=l,4 
DO  12  J=l,4 

ASAVE(K,J)=F(K,J) 

12  CONTINUE 
11  CONTINUE 

CALL  RG(4,4,F,WR,WI,1,YYY,IV1,FV1,IERR) 

WRITE(23,1007)WR(1) ,WR(2) ,WR(3) ,WR(4) 

CALL  DSOMEG ( lEV , WR , WI , OMEGA , CHECK) 

OMEGAO=OMEGA 
DO  5  J=l,4 

T(J,1)=  YYY(J,IEV) 

T(J,2)=-YYY(J,IEV+1) 

5  CONTINUE 

IF  (IEV.EQ.1.0)  GO  TO  13 

IF  (IEV.EQ.2.0)  GO  TO  18 

IF  (IEV.EQ.3.0)  GO  TO  14 

STOP  3004 
14  DO  6  J=l,4 

T(J,3)=YYY(J,1) 

T(J.4)=YYY(J,2) 

6  CONTINUE 
GO  TO  17 

18  DO  19  J=l,4 

T(J,3)=YYY(J,1) 
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T(J,4)=YYY(J,4) 

CONTINUE 
GO  TO  17 
DO  16  J=l,4 

T(J,3)=YYY(J,3) 

T(J,4)=YYY(J,4) 

CONTINUE 

CONTINUE 

NORMALIZATION  OF  THE  CRITICAL  EIGENVECTOR 
CALL  NORMAL (T) 

INVERT  TRANSFORMATION  MATRIX 

DO  2  K=l,4 
DO  3  J=l,4 
TINVCK, J)=0.0 
TSAVECK, J)=T(K,J) 

CONTINUE 

CONTINUE 

CALL  DLUD(4,4,TSAVE,4,TLUD,IVLUD) 

DO  4  J=l,4 

IF  (IVLUD(J) .EQ.O)  STOP  3003 
CONTINUE 

CALL  DILU(4,4,TLUD,IVLUD,SVLUD) 

DO  8  K=l,4 
DO  9  J=l,4 

TINVCK, J)=TLUD(K,J) 

CONTINUE 

CONTINUE 

CHECK  Inv(T)*A*T 

CALL  MULT(TINV,ASAVE,T,A2) 

P1=A2(3,3) 

P2=A2(4,4) 

PEIG1=P1 

PEIG2=P2 

WRITE(25.1008)A2(1,1) ,A2(2,2) ,P1,P2 


C  DEFINITION  OF  Nij 

C 

N11=TINV(1,1) 

N12=TINV(1,2) 

N13=TINV(1,3) 

N14=TINV(1,4) 

N21=TINV(2,1) 

N22=TINV(2,2) 

N23=TINV(2,3) 

N24=TINV(2,4) 

N31=TINV(3,1) 

N32=TINV(3,2) 

N33=TINV(3,3) 

N34=TINV(3,4) 

N41=TINV(4,1) 

M42=TINV(4,2) 

N43=TINV(4,3) 

N44=TINV(4,4) 

C 

C  DEFINITION  OF  Mij 

C 

M11=T(1,1) 

M12=T(1,2) 

M13=T(1,3) 

M14=T(1,4) 

M21=T(2,1) 

M22=T(2,2) 

M23=T(2,3) 

M24=T(2,4) 

M31=T(3,1) 

M32=T(3,2) 

M33=T(3,3) 

M34=T(3,4) 

M41=T(4,1) 

M42=T(4,2) 

M43=T(4,3) 

M44=T(4,4) 

C 

C  DEFINITION  OF  Lij 

C 

L1=YG*((M31**2)+(M21**2)) 

L2=IXY*(M31**2)+IYZ*M31*M21+YG*M21*M11 
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L3=IXY*M31*M21+IYZ*(M21**2)+M=i'YG*Mll*M31+YG*W*{M41**2) 

&  -YB*B*(M41**2) 

L4=YG*  (2 . 0*M31*M32+2 . 0*M21=i'M22) 

L5=2 . 0  * IXY*M3 1 *M32+I YZ* (M3 1 *M22+M32*M2 1 ) +YG* (M2 1 *M1 2+M22*M 1 1 ) 
L6=IXY*  (M31*M22+M32*M21)  +2 . 0*IYZ*M21=t'M22+M*YG*  (Mll=t=M32+M12*M31) 
&  +2.0*(YG*W-YB*B)*M41*M42 

L7=YG* ( (M32**2)+(M22**2) ) 

L8=IXY*(M32**2)+IYZ*M32*M22+YG*M22*M12 

L9=IXY*M32*M22+IYZ* (M22**2) +M*YG*M12*M32+ (YG*W-YB*B) * (M42**2) 

C 

C 

L15=(A11/D)*L1+(A12/D)*L2-(A13/D)*L3 
L16=  (All/D)  *L4+  (A12/D)  *L5-  (A13/D)  =t=L6 
L17=(A11/D)*L7+(A12/D)*L8-(A13/D)*L9 
L25=(A21/D)  =t'Ll+(A22/D)  *L2- (A23/D)  *L3 
L26=(A21/D) *L4+(A22/D) *L5- (A23/D) *L6 
L27=(A21/D)*L7+(A22/D)*L8-(A23/D)*L9 
L35=(A31/D)  *L1+(A32/D)  =t=L2-  (A33/D)  *L3 
L36=  (A31  /D)  *L4+  (A32/D)  *L5-  (A33/D)  ’•'LG 
L37= (A31/D) *L7+ (A32/D) *L8- (A33/D) *L9 
C 

C=CD/(6.0*GAMA) 

L1A=C* (EO* (Mll**3)+3 . 0*E1* (Mll**2) *M21+3 . 0*E2*M11* (M21**2) 

&  +E3* (M21=*=*3) )  +  (l .  0/6 . 0) * (W-B) * (M41*+3) 

L2A=C* (El* (Mll**3)+3 . 0*E2* (Mll**2) *M21+3 . 0*E3*M11* (M21**2) 

&  +E4*(M21**3))+(1.0/6.0)*(XG*W-XB*B)*(M41**3) 

L3A= ( 1 . 0/6 . 0) * (ZG+W-ZB*B) * (M41**3) 

L4A=C* (3 . 0*E0* (Mll**2) *M12+3 . 0*E1* ( (Ml 1**2) *M22+2 . 0*M11*M12*M21) 
&  +3 . 0*E2* (M12* (M21**2)+2 . 0*M21*M22*Mll)+3 . 0*E3* (M21**2) *M22) 

&  +(1 . 0/6 . 0) * (W-B) *3 . 0* (M41**2) *M42 

L5A=C* (3 . 0*E1* (Mll**2) *M12+3 . 0*E2* ( (Mll**2) *M22+2 . 0*M11*M12*M21) 
&  +3 . 0*E3* (M12* (M21**2)+2 . 0*M21*M22*M11) +3 . 0*E4* (M21**2) *M22) 

&  +(1 . 0/6 . 0) * (XG*W-XB*B) *3 . 0* (M41**2) *M42 

L6A= ( 1 . 0/6 . 0) * (ZG*W-ZB*B) *3 . 0* (M41**2) *M42 

L7A=C*(3.0*EO*M11*(M12**2)+3.0*E1*((M12**2)*M21+2.0*M11*M12*M22) 
&  +3 . 0*E2* ( (M22**2) *Mll+2 . 0*M21*M22*M12)+3 . 0*E3*M21* (M22**2) ) 

&  + ( 1 . 0/6 . 0) * (W-B) *3 . 0*M41* (M42**2) 

L8A=C* (3 . 0*E1*M11* (M12**2) +3 . 0*E2* ( (M12**2) *M21+2 . 0*M1 1*M12*M22) 
&  +3 . 0*E3* ( (M22**2) *Mll+2 . 0*M21*M22*M12)+3 . 0*E4*M21* (M22**2) ) 

&  + ( 1 . 0/6 . 0) * (XG*W-XB*B) *3 . 0*M41* (M42**2) 

L9A=(1 . 0/6 .0) * (ZG*W-ZB*B) *3 . 0*M41* (M42**2) 

L10A=C* (EO* (M12**3) +3 . 0*E1* (Mll**2) *M21+3 . 0*E2*M12* (M22**2) 
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&  +E3*  (M22*=t=3) )  +  ( 1 . 0/6 . 0)  *  (W-B)  *  (M42**3) 

L11A=C* (El* (M12**3) +3 . 0*E2* (Mll**2) *M21+3 . 0*E3*M12* (M22**2) 
&  +E4* (M22**3) ) + (1 . 0/6 . 0) * (XG*W-XB*B) * (M42**3) 

L12A= ( 1 . 0/6 . 0) * (ZG*W-ZB*B) * (M42**3) 

C 

Lll=(-All/D)*LlA+(-A12/D)*L2A+(A13/D)*L3A 
L12= (-A11/D) *L4A+ (-A12/D) *L5A+ (A13/D) *L6A 
L13= (-A11/D) *L7A+ (-A12/D) *L8A+ (A13/D) *L9A 
L14=(-All/D)*L10A+(-A12/D)*LllA+(A13/D)*L12A 
C 

L21= (-A21/D) *L1A+ (-A22/D) *L2A+ (A23/D) *L3A 
L22=(-A21/D)*L4A+(-A22/D)*L5A+(A23/D)*L6A 
L23= (-A21/D) *L7A+ (-A22/D) *L8A+ (A23/D) *L9A 
L24=(-A21/D)*L10A+(-A22/D)*LllA+(A23/D)*L12A 
C 

L31= (-A31/D) *L1A+ (-A32/D) *L2A+ (A33/D) *L3A 
L32= (-A31/D) *L4A+ (-A32/D) *L5A+ (A33/D) *L6A 
L33= (-A31/D) *L7A+ (-A32/D) *L8A+ (A33/D) *L9A 
L34=(-A31/D)*L10A+(-A32/D)*LllA+(A33/D)*L12A 
C 

R11=(N11*L11)+(N12*L21)+(N13*L31) 

R12=(N11*L12)+(N12*L22)+(N13*L32) 

R13=(N11*L13)+(N12*L23)+(N13*L33) 

R14= (N11*L14) + (N12*L24) + (N13*L34) 
R21=(N21*L11)+(N22*L21)+(N23*L31) 
R22=(N21*L12)+(N22*L22)+(N23*L32) 
R23=(N21*L13)+(N22*L23)+(N23*L33) 

R24= (N21*L14) + (N22*L24) + (N23*L34) 

C 

P11=(N11*L15)+(N12*L25)+(N13*L35) 

P12=(N11*L16)+(N12*L26)+(N13*L36) 

P13=(N11*L17)+(N12*L27)+(N13*L37) 

P21=(N21*L15)+(N22*L25)+(N23*L35) 

P22=(N21*L16)+(N22*L26)+(N23*L36) 

P23=(N21*L17)+(N22*L27)+(N23*L37) 

C 

C  EVALUATE  DALPHA  AND  DOMEGA 

C 

ZGR  =ZGG(I) 

ZGL  =ZGG(I-1) 

ZGl  =ZGR 
XGl  =XGG(I) 
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c 

CALL  FMATRIX(ZG1,XG1,FF) 

C 

CALL  RG(4,4,FF,WR,WI,0,YYY,IV1,FV1,IERR) 

CALL  DSTABL(DEOS,WR,WI,FREQ) 

ALPHR=DEOS 

OMEGR=FREQ 

C 

ZGl  =ZGL 
XGl  =XGG(I) 

C 

CALL  FMATRIX(ZG1,XG1,FF) 

C 

CALL  RG(4,4,FF,WR,WI,0,YYY,IV1,FV1,IERR) 

CALL  DSTABL (DECS, WR,WI, FREQ) 

ALPHL=DEOS 

OMEGL=FREQ 

C 

DALPHA= (ALPHR-ALPHL) / (ZGR-ZGL) 

DOMEGA= (OMEGR-OMEGL) / (ZGR-ZGL) 

C 

C  EVALUATION  OF  HOPF  BIFURCATION  COEFFICIENTS 

C 

COEFl=  ( 1 . 0/8 . 0)  *  (3 . 0*R1 1+R13+R22+3 . 0=t'R24) 

C0EF2= (1 . 0/8 . 0) * (3 . 0*Rll+R23-R12-3 . 0*R14) 

WRITE  (20 , 200 1 ) XG , ZG , COEFl , DALPHA , OMEGAO , PEIGl , PEIG2 
1  CONTINUE 
C 

STOP 

2001  FORMAT  (7E14.5) 

1007  FORMAT  (4E14.5) 

1008  FORMAT  (4E14.5) 

END 
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